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Abstract 


A numerical method for interpreting. experimental kinetic data 
from multiple-step chemical reaction systems was developed. The 
procedure involves transforming the rate data obtained for each parti- 
cipating chemical component into rate data for each reaction step involved 
in the overall system. The sets of data for each single reaction step 
may then be used to construct an overall kinetic model which is consis- 
tent with the observed chemistry and order of reaction exhibited by each 
of the various components participating in the various reaction steps. 
The procedure may be used with either. concentration-time or differential 
rate data taken under isothermal conditions. 

The method is developed by using. simulated kinetic data contain- 
ing normally distributed random error.. It.is then applied to two sets 
of experimental rate data obtained from the literature. The applicability 
and accuracy of the proposed technique. are believed to be at least compa- 
rable, if not superior, to other. nonlinear estimation methods currently 


available. 
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CHAPTER I 


INTRODUCTION -AND LITERATURE SURVEY 


1.1. Introduction 


Many chemical reactions ocurring in nature or in the laboratory 
do not exhibit a single elementary step between reactant molecules but 
follow a mechanism which involves several such elementary reaction steps 
in various combinations of series and parallel forms. In the laboratory, 
one of the most common methods used to study the behavior of a reaction 
“is to obtain the composition data of the reacting components as a func- 
tion of reaction time. For these experimental data to be useful in 
analysis, future experimental design or simulation, a kinetic model 
describing the apparent chemical reaction steps, and the associated para- 
meters ( here parameter is used for both rate constant and order of 
reaction ) must be determined but consistent with the experimental data. 
In view of the modelling and computational problems encountered in the 
analysis of multiple-step chemical reaction systems: the complexity of 
the chemistry, the various combinations of parallel and consecutive 
multiple reaction steps possible, and the uncertainty in the reaction 
order associated with each reacting ecient. an approach would be 
desirable in which a kinetic model may be determined from experimental 
data, and at the same time, the conventional trial-and-error procedure 
may be eliminated. Some representative multiple-step chemical reaction 


system from recent literature can be described as follow. 
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Tee tks Dehydrogenation of N-propanol 


Wanke [37] has studied the reaction of dehydration and dehydrogena- 
tion of n-propanol on NaOH treated Alundum catalyst at a reaction tempera- 
ture of 463°C and space time ranging from 0 to 150 (hr-gm of cat./gm-moles 
of feed). The reaction scheme is shown in Figure 1. It exhibits a concur- 
rent reaction step of dehydration and dehydrogenation of n-propanol, and 
subsequently from propionaldehyde, two alternative paths, the aldol route 
and ester route, are branched. From the aldol route, through the interme- 
diate product of 3-pentanol, and also from the ester route by hydrolysis, 


diethyl ketone is formed. 


1.1.2. Vapor phase catalytic oxidation of naphthalene 


Ioffe and Sherman [18] have studied the vapor phase catalytic air 
oxidation of naphthalene over V,0. catalyst at reaction temperature 
ranging from 240° to 400°C. The stoichiometry for the exothermic reactions 


which form the various products can be described as below: 


Cj oH + 3/20, ig C1 6°. + H,0 4138 kcal 


naphthalene naphthaquinone 


Ci ots + 9/20, > CoH, 0, oa 2C0, + 2H,0 + 453 kcal 


phthalic anhydride (1-1) 


+ 
Ci ols ar 90, mg C,H,0., + 6C0,, te 3H,0 898 kcal 


maleic anhydride 


ul 
Ci ote te 120, Fa 10C0,, ae 4H,,0 + 1297 kca 


From the reaction system as shown in equation (1-1), and the results 


obtained from their experiments, they have proposed a kinetic model as 


shown in Figure 2. 
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Figure 2 Chemical Model for Oxidation of Naphthalene 


The above kinetic model involves a closed loop series of reaction 
steps, and since the analytic solutions for the differential equations 
describing the kinetics of this system cannot be obtained explicitly a 
priori, more difficulties will be encountered if one tries to use some 


methods involving integrated rate equations. 


1.1.3. Vapor Phase Oxidation of O-Xylene 


The equations written in stoichiometric form for oxidation of 


O-xylene over a SiC-supported V,0. catalyst at reaction temperatures from 


400 to 450°C are 
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Call, + 30, > 3H,0 + CyH,0, + 311 kcal (1-2) 
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Simard, Steger, and Arnott [35] in their study of this system 
have proposed the kinetic model shown in Figure 3. 
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Figure 3 Chemical Model for Oxidation of O-xylene 


The same difficulty as described in section (1.1.3.) will be 
encountered in analyzing the kinetic model shown in Figure 3. 

Before one proceeds to classify various. methods. of estimating 
parameters of a given kinetic model for a multiple-reaction system such 
as those just. described, it is necessary. to. review. briefly the convent- 


ional methods of handling single step reactions, and the accuracy of the 


parameters thus obtained. 


1.2, Methods of evaluating single-step. reaction parameters 


The majority of methods for estimating rate constants and order 
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of reaction for a reaction system involving only a single step, can be 


classified into four categories: 


1.2.1. Method of integrated equation 
1.2.2. Method of differentiation 
1.2.3. Method of half-life times 


1.2.4. Other special techniques 


These methods will now be examined briefly in the following 


sections: 


1.2.1. Method of integrated equation 

A reaction order is assumed with respect. to each reactant and the 
differential form of the rate equation is. then integrated to give a rel- 
ation between the concentration and the corresponding reaction time. 

Take a reaction of the type, 2A > products, which follows a 
second-order rate equation, as an example, using the nomenclature of page 
a2, 

d Cy 


se ce: (235 
de® A 


Upon integration from t=0 to t=t and corresponding C,=€), and 


Ci=Cy> one obtains 


i ea Pe (1-4) 
; kt 1 


The adequency of the integrated equation in fitting the experi- 


mental rate data may be tested graphically by the linear plot of u/Cy 


- o@- ) 7 | 
4 : 


ad aso .qyie s{gute a ylao gatvioval msseye nGitosst & 10% moks2 
teaktogesso ivol osal botttee: 


f 


folisups betetgsiat to bodieM Lt 
norteisoststttb to bodseM o£ 865. 
semti siil-tisd to bodasl .2.Svi 0) 


asupiaross Istosqe ted20 .A.Ssf- 


2i twollot sz at ylisizd benisasxs ad wor | ew shortens seedt 


Ps ta 


:enol3os 
; a . 

f ay 
noijsyps bsitsigeiat to bodjoM .L.8.2 im ’ 


— 
1 


oti? bas tastoses dose o3 Josqest diiw bemuees 8i ssbz0 molzoset A 
-[at s evig oF betaiystnt asdt si soltsups sist sd? to mxo2 istinezei 


-smia noltos91 ae tet bade sdt bas notisiimsones od3 asowsed nok 
s swollot dotdw ,etoubouq +: AS .saqyt edt to nolktoset 8 sist * 9 


P ao, 
sR8q to s1rudslonomon eft goteu -slqmaxs 18 8&5 OLS BUS 9381 zabzo-baeoss 
“oer 
i “re 
i ee 
(€-£) : a ote ba yp. 
id } = F ea ‘ s : 
- ’ , 4a > Spry — 
—_ 
* ‘ . 
brs i: 2 
M8 on 9=,9 gn ‘baoge: sizes bas 7*3 od 03, mort reli = 
. ‘ aa ‘3 F he . 


antasdo 00 «,: 


<4 c 0 i tered 4 ” " 


i 20 ee 


ar. Ne? 3 


4 

<a 

tae age basen 2 
HZ abated sien: hatin 
i Sar % 


ee 
‘ 


: 


against time. Laidler [21] points out that the integration method 
characterizes a reaction according to the way in which the concentration 
of a product or a reactant varies with the time, and this may give rise 
to misleading deductions abogut the way in which the rate varies with 


the concentration of the reactants. 


1.2.2. Method of differentiation 


In this method, values for the left-hand-side of the differential 


form of rate equation, “oe are obtained by graphically or numerically 
dt 


differentiating the concentration-time data. On taking logarithms of 


both sides of equation (1-3) one obtains, 


- log[dc,/dt] = logk + 2logC, (1-5) 


The slope or the intercept obtained by plotting the logarithm of 


rate vs. logarithm of concentration, or by plotting dc,/dt, versus OR 


will directly give the rate constant. 

Using the procedure summarized by Walas [36] of plotting the data 
by hand and then using optical methods to fit functions to the data, 
often unintentionally biased as well as tedious, it is difficult to 
obtain rate data sufficiently precise to establish the reaction order 
correctly. Also, because of the scattering of the data due to experimen- 
tal error, considerable further uncertainty may be introduced by the 
process of numerical differentiation using the conventional 3-point or 


5-point difference formulae. 
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This weakness will be shown in detail in chapter 3, 


1.2.3. Method of half-life times 

The half-life of a reactant; the time required for one-half of a 
given reactant (x=0.5C, 5) to be consumed, is also an useful approach. 
For a second-order reaction, integration of equation (1-3), between the 


limit = = = = i 
imits, t=0, aN Co and t=t, Cy 0.5C os gives 


sel 
Pais [KC 6 (1-6) 


The rate constant can be obtained from the measurements of half- 
life for different initial concentrations. Equivalent techniques are 


used for other orders as well. 


1.2.4. Method of reference curves 


Walas [36] suggests a method of estimating. the: reaction order 
without finding. the rate constant. Here,.a reference plot of ratio of 
the time required to reach any conversion by the. time required to reach 
90% conversion versus per cent conversion is provided. The data are 
plotted on the same scale and on the reference plot, which consists of 
a family of curves corresponding to various orders. The order of 
reaction is estimated by selecting the member of the family of curves 
which fits most satisfactorily. 

Very often, this method becomes impractical, because of the 
requirement of large values of conversion, such. as 90%. Moreover 
since the family of curves on the reference plot: are: limited, only an 


approximate value of reaction order can be obtained. 


- 8 -. ih. 


a 


\ .€ t9e¢qedo ak Its3eb al mworle od {liw aesnisew 


a ettf-tied Yo borteM .€.8.f  —_ 


& to 2Led-en0 xoit bexiups omit sda ;3msj58st & to stil-tisdA eAT 
.dosotgqe futsey os oels ei ,bembenos ad o4 (gad. Oma tastose7t asvig 
eft asswisd .(£-[) sottsups to nolisrgeint ,noljosst sob10-bm0298 & to 
a 
ov 52.08,.9 .t=4 bos . 08.9 ,Oed-,actekl 
oe eS ee ee ee 


Ao 
15 


(8-1 ) 


-tisd to sjnsmemessm 313 mort beatstdo od m&5 .dnsdeno siex sat 
‘g%8 esupiaross tnolsviupa ,anoltistinssnoo istorat tus1stitb toi s tht 


[low e& atebro teito wel 


29viuo so9netstex jo bodiam .+. $ t 
tebro notiosey sf3 gabtemties to bordism s csnguene [oe ] enka ae 
= wes oy % oA, 


to oljs1 to tolq aonsrsis: 6 .sisH .tasienoo 93 x sda gakbat? suodsiw 


oe 


doset o3 batiepet emit sii yd nolasavso5 qos doset of bexiups1 smts oft 
am! > 


- Ss 


sis siab efT ‘.bebivorg al no ketSviieo ‘gm99 19q es notetavaos & 
to watnine dotiw ,telq sanstsI67 oi ‘a0 bas sisoe smme —_ fo beaz0lq 
to tsbao sft .arsbio awotrsv o2 gnibnoqest109 eovaus » to chimed 
geviuo to ylimst ofa to a ond. antiosiee yd betamises: er notsoae = 
Ubresoetetase teom re 

‘ora an seusosd ,Isotsssrqmi semoaed boritem ebay ane “r 


wsvostoM oe as ove nobexswaes to anne: waned to a8 


ud cng S aeciaiaend saneretes oa 60 eeviwo Yo im ” aa 


= (2 


As outlined, even for simple reaction systems the analysis for 
_- determining reaction order and rate constant is very cumbersome by using 
the method of integrated equation when the order of reaction is unknown. 
Although the method of differentiation appears to be more promising, this 
is true only when the accuracy of the rate data obtained by numerical 
differentiation of composition-time data is within experimental error. 
Thus, for complex reaction systems involving various combinations of 
single-step reactions in series and/or parallel form and cyclic or close 
loop combinations abundant problems may be expected. 

During the last decade, a considerable amount of work in the area 
of estimation of rate constants for a given kinetic model for a multiple- 


step reaction system has been done [ 4,16] 
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Before one proceeds to classify the various methods, it is 
desirable to clarify the conditions under which the reaction occurs, 
the mathematical representation of a multiple-step reaction system, 


and finally the method of non-linear least squares. 


1.3. Reaction conditions 


The hypothetical kinetic model which will be presented to 
demonstrate the applicability of this work in chapter 4 and all other 
cases, except where stated otherwise, will be assumed to occur isother- 
mally in a constant volume perfectly mixed reactor. Since the analysis 
of the rate data obtained under these experimental conditions does not 
discriminate between chemical or physical rate-controlling influences, 
to all practical purposes then, the reaction systems are treated as 


though they were homogeneous ones. 


1.4. Mathematical model for a complex reaction system 


For an m-component multiple-step reaction system, the rate of 


each reacting component can be mathematically expressed as: 
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where the terms, veg ys ae denote any general concentration- 

dependent function. 

In compact form, one can rewrite the set of equations by 

1 ij eee : pise 
i 

The rate of reaction for a system which fulfils the reaction 


conditions described in section (1.4), can be expressed as : 


Cees eres = 
ij j ij ae 
Therefore, 
n 
d 
seein yin ey (1-10) 


1.5. Method of linear least squares 


The method of linear least squares which is used as a criterion 
for estimating the rate constants for a multiple-step kinetic model can 
be summarized as follow. Integrating the differential equation (1-10), 
which represents the rate of each reacting component, from Soe to 


t=t > one obtains, 
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If Y, is used to denote predicted values of [C,(t )-C,(t_)], 
iu ty et 


q < represents the experimental values, then the rate constants Ky 
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can be estimated by minimizing D which is the sum of the squares of the 


difference, Oe 


Mathematically D can be written as follow: 
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Since that D is unconstraint it may be minimized by equating 
CBMs = 0, then equation (1-12) becomes 
Q m A 
22 ~ -. J] = - 
2 a [(Y, at 0 (1-13) 
u=1 i=l 
By employing equation (1-13), rate constants a can be estimated 
for a set of experimental data. 


Now, the various methods of estimating rate constants for multi- 


step kinetic models may be classified as follow: 


1.6. Methods of evaluating rate constants for multi-step kinetic models 


The methods of estimating rate constants for multi-step kinetic 


models can be generally classified into six categories as follow: 


1.6.1. Analytic integration of differential equations, followed 
by non-linear least squares regression to obtain values 
of the rate constants. 

1.6.2. Numerical integration of differential equations, followed 
by non-linear least squares regression. 

1.6.3. Differentiation of experimental data followed by linear 
least squares regression. 

1.6.4. Simulation on an analog computer and trial-and-error 
scaling of parameters 

1.6.5. Approximate solutions developed for specific chemical 
reaction systems. 

1.6.6. Elimination of time variable from the differential rate 
equations. 


1.6.7. Quasilinearization. 


| $5. — a ; 


| Sb aut i Rg. 
(Si-£) | dri uk’ “yi ; me 
) | 


galineps yd besimtotm od vent at Iniesdenosay at’ a en a 
a ae 
eeanaed (fff) naieans neds oO = svar 


» ‘ ‘mt 9. 
(€1-£) @ = [ xe\ VEC pl YW}; ©. oe 
| ft otk ut pr 2 fet teu 


hetdamijeae ad aso a ajnssenmoes Pane dl Lf) sebiedon aniwotgns ya 
.sinb I evdomiraqKs to Is2 6 5 
-I3lum toi eJneden0D sis1x got antes ai ebodieg Suckkiee ait, ‘wil = 
swoliot es batitgasts od aie alshoar Shientd qete ; 


efebom slisitd yate-tiium tot siastenod sts1 gntzsuleys ao _pbosieM 0.1 


obisatid qete~iifum tot einsteno» etst gnitsmijes To eborasit off 


swollfo? es eastiogeis. xia osnt betttees ia ylisrensg sd msn aisbod 


o 


bawollot ,anoidsups [stans1siiib to: nokisrgasni ate Debs San 


~~ | 


esulev otgide ot nokeastge1 eszsups tesol rsonki-oom yd 


-ednsdeqes star godt To 


7 Vet 
bewollo2 enol Bups fertnsystitb to nobiexgoact {po brsmutt: 2 veod i 
> 4 _ 
obivanibaics oioas i Jeno! teonkl-n69 wd Asie 
Ate) aut 
sesnil yd bswolfo? atsb tesneattaqxs es notssiinosstitd “Ee a.t \ -— 


‘folsaexget nateupa deaak« < 


t 


torre~bas-Ist13 bas rejuqmos golans 1s Bo notsebuae ae 


ices 9 cia ei ti 


a 


Now, the various methods listed in section (1.3) will be: examined 


briefly as follow: 


1.6.1. Analytic integration followed by non-linear least 
squares regression 

The integration of the differential form, of rate equation per- 
formed in this method is equivalent to the method of integration for 
single-step reaction system described previously: in section (1.2.1.). 
The order of reaction with respect to each reacting component: must be 
assumed a priori. Subsequently the rate constants: are estimated by 
using equation (1-13). 

This method may work successfully, even: with: an approximation, 
with systems involving sets of ordinary. differential rate equations such 
as those obtained for simple reaction mechanisms. But more often it is 


not a suitable method since most realistic models. are not simple ones. 


1.6.2. Numerical integration, followed by non-linear least 


squares regression 
In. this method, the integration. step. is. no. longer limited to 
forms of ordinary differential equations. for: which: analytic: solutions 
are available, consequently, the order. of. reaction can: be either a frac- 
tional number or an integer. 
The partial derivatives of You with. respect ss in: equation 


(1-13) can be expanded as below: 
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Himmelblau et al. [16] and Ball et al. [4] have ignored each 4i-' 


partial derivative term under the assumption of. small contribution by 
of. 
aE Thus, equation (1-13) reduces to 
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epee 


dt-¥, J x [Me jae) = 0 (1-15) 

Himmelblau and his co-workers have assumed the functional relat- 
ionship, ents between time and concentration during integration. Ball 
et al. have applied their procedure in a similar way to a reaction system 
of five components in which each of the reaction steps: follow second- 
order kinetics. To reduce the computations required, they suggested that 
integration of only three differential equations be: carried out for their 
reaction system, the remaining two unknowns can be computed: from two 
independent material balance relations based on the original composition. 

Later, Hunter [17] has shown that for the least squares areyenion 
to give the maximum likelihood estimates of the required parameters from 
multiresponse data, normally distributed data, equal error variances for 
all responses, and absence of error covariance between any two reponses 
were required. 

Eakman [12] claimed that these restrictions cannot be generally 
applied to kinetic data. With the rate constants estimated from equation 


(1-15), he has introduced an additional step involving the hill-climbing 
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method to minimize the determinant criterion. This criterion has been 


derived by Box and Draper= [6], and may be expressed as: 
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The terms which make up this determinant are 


: 
Sag 7 yh, [Oy CED-6, Ce rE: P: ee a aa jet” (1-17) 


on the main diagonal and 
Ss #S - 2 te Ge Ce @ ae ie 2 my k ae 
om yen Oy (tg) =6, (E,) sf, F448 


fe yyat] (1-18) 
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in other positions. Using this criterion, the only condition required 
is normal distributed error around the true value. With a given kinetic 
model, the method estimates the rate constants effectively. But with 
kinetic models in which the order of reaction is unknown, the procedure 
is still cumbersome. 

In the case of complex reaction systems following first-order 
kinetics, Wei and Prater [39] have treated this class of problem rigo- 
rously and in detail. They showed how the linear first-order ordinary 


differential equations (coupling variables, in the case of reversible or 
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cyclic reaction system) can be transformed to an equivalent set of equa- 
tion describing a hypothetical system, each equation of which contains 
only one dependent variable. The equivalent set may then be solved by 
standard matrix algebraic methods to give the usual exponential type of 
integral equation as the solution. With this matrix transformation, they 
successfully reduce problems of the class in section (1.6.2.) to the ones 
in section (1.6.1.). 

In this method, to perform the transformation, the square matrix 
K of unknown rate constants must be used and consequently, an iterative 
calculation procedure is required. In spite of its many elegant features, 
the Wei and Prater method is limited to first-order kinetic models. 

When multiple-step reactions can be approximated by a first-order 
model they may be handled by their method, but very often this approxi- 
mation is satisfied only for a limited region of the experimental condi- 
tions. Ames [1], Aris [2], and Wei [38] have independently shown how 
a non-linear reaction system can be reduced to a linear one, and thereby 
the system in the linear form can be handled by Wei and Prater method. 
But the mathematics of transformation is rather complicated and the met- 
hod is still limited to system exhibiting simple non-linear kinetic 


behavior. 


1.63. Direct differentiation of composition-time data 


Levenspiel [24] describes the common graphical procedure showing 
how the reaction rates for a set of concentration-time data, can be 
obtained. The procedure involves smoothing the experimental data points 


first by drawing smooth line through them and then tangents were drawn 
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at several points on the smooth curve corresponding to different time 
intervals. Subsequently, linear least squares regression was applied to 
the rate and corresponding concentration data. If the order of reaction 
is other than first order, a logarithmic transformation procedure, as 
described previously in method of differentiation for simple reaction, 
was applied. In this way, rate constant and order of reaction can be 
estimated at the same time. 

Because the method of smoothing of the experimental data which 
was used and the graphical method of differentiation ais generally 
influenced by the experimenter's biases, the results can occasionally 
be misleading or ambiguous. Its outstanding feature is of course that 
rate constant and order of reaction can be obtained simultaneously. 

Bak [ 3] recommends a procedure for smoothing the experimental 
data, measured at equidistant values of time, using third order ortho- 
gonal polynomials to approximate the data by minimizing the sum of the 
squares of the differences between the predict values and the experimen- 
tal data values. He then proposed to replace the experimental data 
points by the corresponding values of the approximated polynomials. The 
reaction rate is obtained by differentiating the polynomials. Subsequen- 
tly, rate constants can be approximated by applying linear least squares 
to. fit the data. 

Lindsay [25] observed that one of the reasons why the reaction 


rate,dC,, cannot be obtained precisely is due to the variation in magni- 
dt 


tude of dCy over the course of an experiment, i.e., the proportional 
fat 
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error:- the ratio of the error in 


Tent to Teac In order to minimize 


the error, he proposed a procedure for finding an optimum time schedule 
for making composition Rosaurementst By transforming the time variable 
to a dimensionless variable, so that the changes in Cy between one 
sampling time and the next should be equal throughout the course of the 
experiment, the experimental error was minimized. In principle then, 


accurate rate data may be obtained from composition-time data. 


1.6.4. Analog computer, trial-and-error approach 


The rate constants are estimated here by simulation of the 
mathematical model for the reaction system on an analog computer. The 
procedures can be outlined as follow: 

(1) Construction of a block diagram for connecting the operational 
elements corresponding to the differential rate equations which described 
the assumed reaction mechanism. 

(2) Selection of the scales for representation of the dependent 
variable and time. 

(3) Calculation of parameters of the model from the coefficients 
of the original equations and the selected values of the scale. 

The criterion generally used, is to select the mathematical model 
whose characteristics compare most favorably with the experimental curves 
for the composition of the reaction mixture as a function of time. The 
curves simulated by the analog computer provide a visual evaluation of 
the alternative models. An example of this, is given by Katarov and 
Lutsenko [19]. As ae criterion employed in this method, is to compare 


the shape, the major deficiency of the method lies in discriminating 


tmiotens 1? Ya .ednenmeivessm nolLiteoqmos ‘gif “x03 *~ 


(25 
t : f t ; “_ os > t . = 
goltbaogmos moit bontside sd Yam BIRd S387 sINTLO Ss 
i” 


mw on A 


soulayv betjosisa oda bos anoltisups fe 


.smtt to noldonv? s es Siudxim noksoaex odd to nokzteoqmos eda x 0% 


ee 
16 maiouls vo Isuztv 8 shivetg xrasugMoa golans ori — Sepatiias o rua 


= Je 
Am 
f La 
4.8 
+ ° = 


> od =————— jt. rorss afi} to obdas ota ayn 


Pooky TOY orrube jo 6 its 
garbok? tot srubsso01q & bseogorg on r0738 


gadsiio sd) tata of ,oldstrev isstackeuaiall 7 


orii Isupe od biuode dxen sd3 base om 
a. 
sesimiatm eaw yorrs letoemiteqxs ad2 ,Jneatseque 


- 
, - . * # 
so%yve-bos<[etto . trqmoo solenmk -#.d.1 
a _ 4 — — _ ee te I = 
sd sips bsisniies sta eineJenos eiex ont i 


>’ BF 


nottosex sia tol Jebom Isolinme at 
men’ 


gje£7 intineausiitdb sria G3 gninnogestsz05 alana i. 


.usinadosm noktoset beaues oid 


9297%q! 2 asigoa sii to soksosle® (s) - 
Sait bone of 
Isbom adit to eipvemeitaq to aotsalunled (€) 


isa of al ,bsayw vilaransg wonton oat 


iol 
siw edszovst Jeom exaqmoo eoise trospanado 99¢ rhe 
mn a7 a 


‘s , oy _— 


ae 


: «ahd? Io 
ae > 


erate 
between various "best fits" for alternative reaction mechanisms. 


16-5. Approximate solutions for a specific reaction system 
Because many addition, substitution and polymerization reactions 
take the following irreversible consecutive, competitive form: 
Ki 
CVE nL a a. ee er 
1 
ae 
Ter B f Py “h e oe * 6 (1-19) 
+B +3 | ee 


Po 3 e ° ° 


and because the general solution to the system with any initial com- 
position and any set of constants is still not available, Friedman and 
White [14] proposed a procedure of transforming the differential rate 
equations which describe the reaction system, into dimensionless form 
so that the results would be applicable to a wider variety of systems. 
They introduced an empirical parameter b which is a function only of 


K, (=Ki ) to characterize the relationship between B and te by the 


ky 


below equation, (1-20) 


8 = _e > tt (1-20) 


: 
Then they derived the general approximate solution to the system 
from the ‘differential rate equations which represent the reaction system 


in equation (1-21) 


Peactan 
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In equation (1-22), an arbitrary initial value of 8 is related to 


the empérical parameter b which is a function only of K (aK ye Ror 

i. iy 

Soa or. greater than 0.1, the values of. parameter b can be calculated 
from a second-degree multinomial in “s and EES The nature of the 
dependence of parameter b by the approximate solution is found by fitting 
the predicted values calculated by equation (1-21) and the experimental 
data values through the method of least squares. 


This approximate solution method handles well the reaction system 


in equation (1-19) but lacks generalization. 


1.6.6. Elimination of time variable from the differential 
equations 
Benson [5], proposed.a procedure for eliminating the 
time variable. by. dividing each of the differential. equations. which 
describe the. reaction system by one.of.these equations. For parallel 


reaction steps such as, 


Figure 4. A Parallel Reaction Steps Chemical Model 


suppose that the order of reaction for step.1 is nD and for step 2 is 


Nj: Then, the rate equations for the two steps may be represented by, 
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Benson has shown that dividing equation (1-23) by equation (1-24) 


to eliminate the variable yields 


PGperabarl zt 
Bolte cl 


37") >) (1-25) 


rf 144° Ni? then direct integration is possible, yielding 


k 


Cc, =—C,+9 (1-26) 
2 ki i 


The integration constant, $, is usually zero at zero time. If C, 


and C3 are produced from identical reactants in which the overall kinetics 

are of equal order, they will be formed in a constant ratio at all times. 
Mezaki [27] points out that because of the experimental error 

encountered in measuring concentrations, the rate constant calculated by — 


using this method may be erroneous. This is especially true when the 


reference concentration C, of the component cannot be obtained accurately 


(in equation (1-25)); in this case, the rate constant estimated will 
deviate very much from the true one. The other obvious disadvantage of 
using this method is that the rate constants can be obtained only 


relative to one another. 
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1.7. Mechanistic modelling using a set of experimental rate data 

In the last section, different methods of estimating rate constants 
for multi-step kinetic models have been discussed. Not much progress 
has been accomplished in the last decade on the problems of mechanistic 
modelling, in the strict sense from a set of rate data. In this section 
some of the most recent works will be examined briefly. 

Peterson [29] has outlined a procedure for the classical two-step, 
consecutive reaction system with an unassigned order of reaction, M1. 


to the reactant, A, and a first order exponent for the second reaction 


uy 
step, 
F Sel Koo . 
1 2 3 
Ld 
Figure 5 Consecutive reaction system 
The differential rate equations are written, using the same 
nomenclature, 
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By using non-linear estimation in conjunction with analytic 
integration, i.e. the method mentioned in section (1.6.1.), the rate 
constants and order of reaction can be approximated. He has shown how to set 


up an intermediate mechanism as a stepping-stone towards building a more 
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adequate. mechanistic rate equation. for. the complex reaction system. The 
standard deviation of the fit and: physical reality. of the results (ex- 
pressed as. positive values for the rate. constant) are the two criteria 
used to discriminate between the. various: plausible models. 

Fux and Ioffe [15] in their. study. of.mechanistic modelling claimed 
that because of the unknown nature. of experimental error as well as the 
unknown. values of rate constants for: a real. reaction system, it is more 
reasonable to carry out the initial analysis, using a set of simulated 
data, which characterized the actual. kinetic experimental points stat- 
istically. At first, true data were. generated from a set of differential 
equations. which describe an assumed. monomolecular kinetic model in which 


all the reaction steps follow first-order kinetics, 


Figure 6 Consecutive and Parallel Reaction System 


with. initial condition at t=0 C,=1, C,=C,=C,=0, and with appropriate 
values for the rate constants. Then normally distributed random error 
with. standard deviations, s=0.02; 0.05; 0.1; 0.15, respectively, were 
introduced into the true values. .To. demonstrate that the true kinetic 
model can be obtained from the experimental data using their procedure, 


they have assumed several alternative reaction schemes for their simu- 
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lated "experimental data". The rate constants.of a kinetic model were 
estimated first by integrating the differential. equations which represent 
the. model. and. then: minimizing the: sum. of. the. squares of the difference 
between the integrated values and the: simulated "experimental" ones. 

The adequancy of the assumed reaction scheme in fitting the data was 
checked by using only positive values. of the rate constants as the cri- 
terion. If negative rate constants were obtained, the estimation of 
different values for the negative rate constant was repeated by using 


simplex method. They have shown a monomolecular reaction scheme of the 


type, 


Figure 7 Monomolecular reaction system 


can be reduced. to the true kinetic model of. Figure 6 from the simulated 
"experimental" data. They also suggested: that if an independent exper- 
imental study of a particular reaction step. were performed, improved and 
more consistent results can be obtained by. calculating the rest-of the 


rate constants in terms of the independently measured value. 
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CHAPTER II 


ANALYSIS OF MULTI-STEP CHEMICAL REACTION SYSTEM 


2... Introduction 

As pointed out in chapter 1 due to the complexity of the chemistry, 
the various combinations of parallel and consecutive multiple reaction 
steps possible, and the uncertainty in the reaction order associated with 
each reacting component, an approach would be desirable here in which a 
kinetic model may be determined from experimental data so as to reduce 
the modelling and computational problems encountered in the analysis of 
multiple-step chemical reaction systems. A semi-empirical approach 
represents an improvement over strictly empirical methods for developing 
the model, jf one recognizes that kinetic studies of much simpler 
reaction systems than those cited earlier have not yet departed from a 
largely empirical correlation point of view. 

Numerical differentiation of composition-time kinetic data has 
attractive features in analyzing complex reaction systems. With the 
customary "concentration-time" experimental data, Dalla Lana [10] has 
suggested that the individual steps in multiple-step chemical reaction 
system be modelled. This would reduce the degree of empiricism associated 
with procedures requiring direct modelling of the entire system. The 
modelling procedure can be developed by considering the general form 


of rate equation shown below, 
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in which the contributions to the overall rate of change of one compo- 
nent from each reaction step are summed. 

The individual terms in the right-hand-side of equation (2-1), 
"43? can each be approximated by a power type. of rate equation. This 
approach is commonly followed in the analysis of simple one-or two-step 


reaction systems, and gives for each term, 
i pete tata Oh, (2-2) 


for a monomolecular reaction step, Ay > product 

Local values of the left-hand-side of equation (2-1) are known 
from numerical differentiation of: concentration-time experimental data. 
When the number of independent ordinary differential equations ¢ @fethe 
type, equation (2-1), is equal to the number of contributing steps, then 
local values for each reaction step, 45? can also be expressed in terms 


of the known values of oe By this approach, the analysis decomposes 


dt 
overall rate data into rate data for individual reaction steps. The use 
of simple algebraic computations required in this approach to modelling 
is more attractive computationally than the use of variations in the 
integral method when the set of differential. equations has to be solved 
simultaneously to check the assumed model. 

The remainder of this section. will. provide the groundwork for 
the subsequent applications of the proposed techniques. Before trying 
to demonstrate the applicabability of this approach, a brief survey of 
numerical differentiation methods, and a method of obtaining accurate 


first derivatives in the situation when a small amount of data are 


available will be presented. Following this, time--composition data 
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generated from a numerical solution to a hypothetical non-linear kinetic 
model, plus normally distributed random error will be used as a first 
demonstration to examine the applicability and the accuracy of the 
proposed approach. To complete this study, the method is applied to 
several real reaction systems, the data for which are reported in the 
literature. 
2.2. Classification and analysis of complex reaction systems 

The system of algebraic equations, which must be solved for 
various kinetic models can be classified into one of three cases, 


depending upon the "chemical model” involved: 


2.2.1. Open network reaction systems (intermediate or final 
product components do not form from two or more converging 
reaction step): 


Open network examples, 


A; — Az A, 5 w 
A 


2.2.2. Closed (or loop) reaction system (intermediate or final 
product components may form from two or more converging 
reaction steps): 


Examples of closed (or loop) reaction systems, 
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2.2.2.1. Reversible reaction system (subclassifcation 
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In the analysis immediately following this, the reaction rate 
dc. 
component, i(=r;), is assumed to be known from the differentiation 


dt 
of the experimental composition-time data. 


oo 2 se Open network reaction. systems 


Suppose a reaction system follows: the: chemical model shown 


below: 
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Figure 8 An Open Network Chemical Model 


dC. 
The reaction rate components, 1, (i=l,. . .,5), all measured at 


dt 
a specific reaction time are each equivalent to the algebraic sum of the 


rates for each contributing reaction step, ae Mathematically, the set 
of sore can be expressed in the following form, where a is now used as 


an abbreviated notation for ane 


dt 
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Since all r, are available from the differentiated experimental 
data and since only four unknown reaction steps. occur on the model 
shown in Figure 8, for a common "time", local values of rates for each 
reaction step, Ts? can be obtained simply by solving the set of equations 
(2-3). This procedure is used to generate a set of r j at a number of 
different times, in effect, providing: a set of. rate data for each inde- 


pendent reaction step. The chemical.model of Figure 8 may be completely 


defined upon estimating all Te a aree by regression techniques. 


2.2.2. Closed (or loop) reaction systems 


The rates of overall reaction for a reaction system with the 
chemical model, 
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Figure 9 Chemical Model for a Closed Reaction System 


can be described mathematically, at a common time, as shown below: 


r,=r,,-° 2 (2-4) 
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Since the determinant of the coefficient matrix, M, is zero the 
matrix is singular and the inverse matrix, (m7, does not exist. The 
Ts column cannot be solved explicitly. Farkas [9] has solved asimilar 
loop reaction system in which all the reaction steps follow first-order 
kinetics, as described in chapter 5, by assuming values for the relative | 
rate constants for the branching reaction rate steps. Of course, the 
values of rate constants obtained in this way will not be absolute values 
and further, similar values for the reaction order of both branching 
steps are required. Experimentally, this problem can be overcome by 
performing an independent kinetic study for one of the intermediate 
components under similar reaction conditions. 


In this work, for kinetic models of loop reaction system it is 


proposed to estimate the parameters (the rate constant and order of 
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reaction) of two reaction steps by non-linear estimation. This can be 
visualized by considering the component, Ao» in Figure 9. The first two 


equation of (2-4) yield, 
r,tr,.=-r - fr (2-7) 


Then, using power function approximation to the two “ig one 
obtains, 


Noo n 
ap bas ety 
Boia 4 kK 2% Koh 2 


24 (2-8) 


Since (r,+r,) and C, are known values, the parameters, k 


2 ye ty 


Koy? and n can be obtained by non-linear estimation. Once the r 


24’ Ze 


and yA parameters are known, the remainder of the reaction steps, a? 
in equation (2-4) can be evaluated by algebraic substitution. The 
components with higher concentration values are usually chosen for the 
non-linear estimation. The reason for this preference will be demon- 


strated in section (4.4.1.). 


2.2.2.1. Reversible reaction system 


Figure 10 Reversible reaction system 


Rate balances around each reacting component are given by: 


(2-9) 
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In matrix form, the equations become 


r 2 att 
1 11 
* (2-10) 
2 02 1-1 
Y=er .M i 
sine" (2-11) 


It is obvious that the determinant: of. the. coefficient matrix of 
equation (2-10) is zero, therefore ae and a cannot be evaluated expli- 
citly. Both rates can be determined only when there is an independent 
means of specifying the value of one parameter or the ratio of parameters. 
The thermodynamic equilibrium constant,. preferably experimentally deter- 
mined, might provide such imformation. In the absence of knowdédeeaabuit 
the equilibrium constant, the rate step can be determined by applying 


the same procedure as described: in: section 2.2.2., i.e. estimating k_, 


it 


k >? and n,, from either one: of: the. equations in (2-9) by non-linear 


mit? oc 2 


estimation. 


22 


The validity of this procedure for sections (2.2.2.) and (2.2.2.1.) 
will be examined by applying this: approach to.some actual experimental 


reaction data. These will be presented: in: the final chapter. 


2,2.3..Combination of reversible: and loop. reaction system 


For a reaction system of the type, 


op. 
ae 
A 3 Figure 11 Reversible and Loop Reaction System 
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the chemical model may be determined by first applying the non-linear 
estimation procedure to the reversible step as outlined in section 
(2.2.2.1.), and subsequently analyzing the remaining reaction steps of 
the system as described in section (2.2.2.). 

The method of analysis of the data described in section (2.2.) 
has little value if the accuracy of numerical differentiation of 
experimental data cannot be assured. A brief survey of the available 
methods for numerical differentiation, the importance of smoothing the 
experimental data before differentiation, and a method of obtaining 
correct first derivatives from a limited set of experimental data will 


be presented in the following chapter. 
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CHAPTER. III 


DATA SMOOTHING AND NUMERICAL. DIFFERENTIATION 


3.1. Introduction 

The measurement of experimental: data is. always subject to errors. 
These errors, may be sufficiently large sometimes that the algebraic 
sign of the first derivatives, obtained: by differentiating the data 
using a simple difference technique, is incorrect. This situation can 


be visualized by considering the example of Figure 12 taken from Ralston 


[31] 


© Experimental data points 


_ True function 


yi wan Exact polynomials 
vis approximation 


Z —w— Linear least squares 
/ approximation 


—-—Exact linear approximation 


Figure 12 Least squares and exact approximation 
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The five empirical points employed by Ralston possess, in fact, 
a linear relationship between X and Y. Since the errors in the data 
are quite significant, linear approximations using any two points will 
lead to very poor approximations to the true function. The sign of the 
first derivatives, obtained by linear approximation, will be different 
from that of the true function. Although the fourth order polynomials 
approximation (exact fit) does not deviate very much, the first deriva- 
tives are very different from the comparable derivatives for the true 
function. In the case shown, the linear least squares approximation lies 
close to the true function with a similar slope. 

Thus, it is essential to smooth the experimental data before it 
is differentiated. Walas [36] and [3] also stressed the important role 
played by data smoothing in numerical differentiation. With the advance- 
ment of digital computer methods, it is no longer justified to smooth 
the data by hand. Furthermore, though the least squares method works 
well in many cases, it does not always produce a satisfactory representa- 
tion of the data. Ralston [32] claimed that a non-linear function can be 
approximated well by a straight line over subranges of the data. This 
is true only in the presence of a large number of data points so that the 
linearity within the subranges is well assured. Lanczos [22] suggested 
a procedure of "movable strip" or moving arc , i.e. evaluating the 
derivative at each center point and moving stepwise down from the top 
of the data. In this Ce eae the first derivative can only be obtained 
accurately at the center point of the moving arc. Suppose a five-point 
moving arc is used: then there will be a loss of four data points (the 


first two points and the last two points) for the first derivatives at 
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these points cannot be evaluated accurately with this approach. The method, 
as a result, may mot be suitable for correlating kinetic data (time- 
composition) of complex reaction systems because the data are generally 
obtained in a limited amount. Thus, the problem becomes how to select 
a function, generally nonlinear, that will approximate the limited amount 
of complex reaction kinetic data (say seven or eight points). 

Many types of functions are available for fitting purposes. In 
the following section, a brief review of some methods employing the most 


commonly used polynomial functions, will be presented. 
3.2. Polynomial Approximations 


For polynomial approximations, the following three general types 


can be classified: 


3.2.1. Interpolation approximations 
3.2.1.1. Newton's interpolation formulae with equal spacing [8] 
3.2.1.2. Lagrangian functions [8] 

3.2.2. Least squares polynomial approximations [23] 

3.2.3. Minimization of the maximum error polynomial approximations 


[32] 


3.2.1.1. Newton's Interpolation Formulae 
For a set of equally spaced h, ee Xo Xd points, 


using Newton's forward difference formula, the function f(x) can be 


approximated by p(x) as follow [8]: 
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£(x) = p(x) = £ + (AE + (A enced Beal OES e (3-1) 


where s = 


ABT Eyew Ege ee £Q Py ove f= PL) 


The crucial weakness of this method lies in the inverse proportion- 
ality of round-off error and truncation error to the equal spaced interval 
h, i.e. by decreasing the value of h so as to reduce the truncation error, 
the error caused by round-off will magnify. Although an optimum choice 
of h may sometimes be made when the two error-bounds are equally good, 
the dominant error source lies on the input errors themselves. From the 
numerical differentiation point-of-view, as pointed out in section (3.1.), 


exact polynomial approximation may not be suitable. 


3.2.1.2. Lagrangian Functions 


This method is useful when the data are not equally spaced. 
Using Lagrange's interpolation formula, the function f(x) can be 


approximated by, 


n 
f(x) ~ p(x) = 52) L, ()y, (3-2) 
n 
Tl (x-x.) 
j= 
3#i 
where, L, (x) =a <r (3-3) 
n 
- HasBy? 
j#i 


are called Lagrangian interpolation polynomials 
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The disadvantages of using this approximation are: 

(1) the coefficients of the polynomial L, (x) are dependent on 
each other, i.e. the addition of one more term to the polynomial will 
require a complete recalculation. 

(2) values of L(x) are obtained by taking the difference between 


data points very close to one another. 


3.2.2. Least squares Polynomial approximations [ 23] 
Suppose a function, f(x,), with discrete data points, i=l,2,.. 
. « » MN, Can be represented by equation (3-4) 
m 
ou apy, + e, (3=4) 


f(x, ) sn 


The principle of this approximation is to minimize 


z 
e 
i=l 


so as to yield best estimate values: for a 


3.2.2.1. Simple Power Functions 


Tho thiss casecaipower series is substituted for p (x,) so that 
4 


equation (3-4) becomes 
r : (3-5) 
£(x,) Pykg RP + By 


where B, are the coefficients: of the power series 


Lapidus [23] has shown that by minimization of the sum of the 


tot Nepean FS) 
squares of the error terms, i.e. if s_ ia -. > 0 (3-6) 


Ones 


ee 


ne B) Ata ; 


“26%8 colismixor 'qas ayes gnteu to pio 


—_ 


no Sasbasqeb ste (x ) Latmonytog ind to etnotolitsos ott wo 


o3 miss. 210m sao La noittbbs oft .a.k vibe, 


aa 


dhby Inimonylog si3 
10k eluslsost ato iqioo 6s 

i 

_ is 


anisei vd beatstdo sts (x), fai: eoulav oe ’ 


meswisd sonsieltib era g 
ie 


.reftjons 9eno.o2 cues roe exato 


4 


a 
BE 


(es ] anotssmtxorggs isimon ylot eorsups eset eS 7 G4 
jntog sisb ste1seib dtiw .(,x)2 <0 f39nut nell / 
: “ 


(8-£) notsaups vd beimsestqex. od MED fh ¢ 
od a 


x . a er 
(d=€) a pr (#48 Se (3 7 


osimiatm ot ef noltsmixoxqqs atdi to si 


atec 


ag, coet £9 (3-7) 


. The following set of simultaneous equations written in matrix 


form results, 


To Ts 
X XB = xy (3-8) 
YW wx. eS é Ix. 
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B is the vector of the unknown coefficients, ve 
The power coefficients, B, are usually determined by Gauss 


elimination. 


3.2.2.2. Orthogonal Polynomials 


Orthogonal polynomials sometimes are a useful tool for approxima- 


ting a function because of the following features: 
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(1) the solution of the linear equation is guaranteed 
(2) the coefficients of the polynomials are independent of each 
other, i.e. increasing the power of polynomials does not 
change the previously calculated coefficients [40]. 


Using an orthogonal polynomial Q (x,) to approximate £(x,) for 


j 


discrete data points (i=1,2, . .., n), the use of a least squares 


criterion gives 


m 
£(x,) = i (x,) + e; (3-11) 


YQ 
j=0 ? 
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In the same manner as before, the least-squares minimization 


leads to 


WV 
i=) 


LOPS 
min S =. (3-12) 


When numerically differentiating relatively smooth equal-spaced 
data, gram orthogonal polynomials proved to be quite effective and 
convenient [40]. The coefficients of the gram orthogonal polynomials for 
evenly spaced data points are previously calculated and the need to solve 
a set of linear equations is eliminated. When the data points are smooth 
and exact, a second-degree polynomial or parabolic approximation is 
found to be adequate. A method of testing the quality of the fit and a 
formula for calculating the standard error of derivatives attributed to 
the experimental error have eo developed [40]. In practice, irregular 


data points give fluctuating and inconsistent results. 


3.2.3. Minimization of the Maximum Error Approximation [32] 
As the title of this sectjon expresses, the basic idea involved 
in this method is to minimize the maximum error that arises due to the 


use of an approximation. The error tern, ess for discrete data, Xi Vy> 
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where i=1,2,3, . . . ., n, can be expressed as shown below, 


en i) a4 (3-13) 


If we assume £(x,) is the function that approximates the data 
points and if E is the largest error among ej» then the object of this 
method is to determine f(x, ) for which E is smallest. The minimum 
maximum approximation is also called a Chebyshev approximation. 

As for what type of function to be chosen for minimum maximum 
error approximation, Ralston [32] presents a full account of the more 
accurate results obtained with a rational function which is made up of 
the ratio of two polynomials rather than the original polynomials. He 
has also developed a technique for generating a rational approximation 
which, among all rational approximations with the same degree polynomials 
in the numerator and in the denominator, gives the minimum maximum error. 

The method of smoothing cited in section (3.2.1.) are unsuitable 
for numerical differentiation. Although least squares method (section 
(3.2.2.)) smooth the data reasonably well,the first derivatives evaluated 
from the resulting least squares approximation sometimes give erroneous 
values (this can be seen in section (4.3.)). 

In view of the following problem: the non-linearity of time- 
composition relationship and the availability of only-limited amount 
of data, it will be advisable to fit a linear or non-linear function 
which essentially has a shape similar to the set of time-composition 
data. Subsequently, the derivatives are obtained by analytical or 
numerical differentiation of the resulting function. 

In relation to this, Mezaki and Kittrell [28] have shown the 


advantages of fitting a non-linear function to a set of conversion-space 
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time data which have non-linear relationship. They fit 
W 
x = Atanh[B() | (3-14) 


to the data, and analytically differentiate the function at the origin. 
The parameters, A and B, define the function which was fitted by non- 
linear least squares. They found that by including the high conversion 
as well as the low conversion non-linear data, they could discriminate 
between two rival models. On the orther hand, fhe discrimination cannot 
be done if the linearity of the data was assumed. 

In this work, the parameters of the non-linear function fitted 
to the data are approximated by minimizing the residual sum of squares 
using Rosenbrock's Hill-Climbing Method. [34]. 

A brief review of the Hill-Climbing method and some of its modi- 


fication will be described in the following paragraph. 


3.3. Optimizing by Hill-Climbing Method [34] 


The optimization of a function arises in many ways in engineering 


application. Some of the common methods can be summarized as follow: 


3.3.1. Steepest Ascent or Descent Method 


In this method [34], the partial derivatives of a function f(x); 


Xoo -, X_) are first found by using a set of p experiments, i.e. 
si 2 x is first varied by an amount of Os» then = can be approximated by, 
hi 


f(a, + a> a By eee had hire f(a, a> re ey a) 


Ox) x=a a, 


(3-15) 


Subsequently, the remaining derivatives are estimated in the same 


way. Then the direction of steepest ascent is found by calculating the 


-! 


pee 
322 yod .qidenotssiss resatl-non aved | ote w stat > smis 
2 4 etic ) a ae 


@. 2a 
(b5-€) } (Gp aldaetA =< : 


-Sigite ef? Je soldony? aris stetiasralithb Yilesliayians boa ,stab aed aa 
gd i ’ doidw nokitonaui sdt antish .& bas A .ey 2s TEq aes 
ieva02 fgis sit gaibutont yd tsd? bavol yesd? -eetaups tenet 38965 


stanimitoelb biuon vada ,.steb reaankl-non nokexravnaeo wol sat aa Lisw 9 
7 { 


* - a 
itmizsetb of3 ,bosd ryeh320 and nO .elebom feviy ows csswis¢ 


ian) 


-bomusas esw sish sti to vituasnil oft if sm: 
bS33i2 soitonut tesati-non sdJ to etsismsxsq of3_,drow etda Al i 


. . h ted 
sersupe io mue Leublesy ofi gatstmtnim yd bstemtxorqae sxe a 


f «ar 


LSE] bodjaomM ant ImE LO-LLITA e'd ,ordapeodk | 


atte 
bos bodsem gatd@rlo-Lith edt to wetves Jsixd A 3 


ee * 


. « , . a) - 
-Agsigsrsq gotwodlot oft at beditoesb. sd Eite aokjeait 


LAE} bodseM gsebkdimt HO~TLiH yd gnisimksqO .£,€ ; 


sotjoaw? & 10 nottestmisqo sAT - 


‘Wollot 25 hesivamme sd nso abodjem nommos sid Yo smoe -nolisotiqys 
SorteM tna aaa 10: d095RA tas 1gae32 -L.€.€ a - 


. -*)2 cokions? 5 io esytisviasb Isksteq sid « dé] Godden sd at 


i 
insmiasqxs q to 35@ & gntey yd hauot Jett orn is a 4 © 0@ 
> a» 


. 
76- 


a 
«vd hetamixoxuqqs 6d nso sr aan r lo Inpyome a8 XG bebisyv dextt at? 


Sc 


unit vector: L with components 
Te af/ax; 


o£ 2,1/2 (3-16) 
C2 Fax) 71 


When p=2, the procedure can be illustrated in Figure 13. 


x 


1 


Figure 13 Procedure of Steepest Ascent Method 


then the highest point b, on the line- estimated from the initial 


point a, in the direction of L, is found. - The whole process is repeated 


again by treating b as an initial point. The search is carried out until 


the summit, d, is reached. 


3.3.2. Rosenbrock's Method. of Hill-Climbing 


Rosenbrock [34] modified an accelerated Hill-Climbing method by 


introducing a variable scaling factor. After each new function evaluat- 
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ion, the values of the function, f, so obtained is compared with the 
best previous value. If the operation was successful in increasing the 
function value, the value of the step length, e, is trebled in the 
favorable coordinate direction. If it was unsuccessful, e is multiplied 
by -+ . The effect of this modification is that if a change parallel 
to u, orthogonal unit vector, (i being a component of a set of p 
orthogonal unit vectors) was successful, then the next trial is made in 
this direction it will be in the same sense and three times as long. If 
the trial was unsuccessful, then the step will be made in the opposite 
sense and will have half the length. So the operation of searching the 
optimum point in fact tries to explore on each line, instead of taking 
a fixed step in each direction. 

Although the progress of Rosenbrock's method may be slow when 
there are narrow ridges, because of its generality it is relatively 
difficult to defeat. In this work Rosenbrock's method of Hill-Climbing, 


which does not require analytic first derivative, was used. 


3.4. Minimizing the Residual of Sum of Squares 


The choice of an approximating function has to be made on the 
basis of experience and intuition [33]. The conventional resolving 
power of polynomials or trigonometric sums are not practical for 
numerical differentiation. Even though they might represent the 
experimental data faithfully. 

Some common functions, Y=f(x), that generally have similar shapes 
to the ah eaehineds toned kinetic data, were drawn by an IBM 1627 plotter 


using the appropriate values of parameters of the functions. These 
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curves are as shown in Appendices (1) and (2). Their use helps to 
choose an appropriate function that will best serve the particular case. 
Sometimes, the accuracy of approximation can be improved simply by 
adding another identical term to the previous function, as illustrated 
in the case of 2 and 3 in Appendix (1). Generally, if the fitted 
function behaves reasonably well compared to the experimental data, 
three to four parameters are accurate enough to determine the approximate 
function. The use of too many parameters approaches the case of high 
power polynomials and at the same time, computing time increases 
immensely. Computing time here in general is a function of the total 
number of data points, the number of the parameters in the fitted 
function and the values of the initial guess introduced. The computing 
time for a function of 4 parameters to approximate a set of data (say 15 
points) with initial guess of unity in an IBM 360/67 will vary from 2 to 
4 minutes. 

At the end of the searching operation for a fitted function, the 
variance of the fit and the deviation between the observed experimental 
data and the calculated values from the fin cetod are recorded. The 


variance is defined as: 


N 2 
a = 7 Casopaer. = St,caten.”” 


N=P (3-17) 


The variance was primarily used as a measure of the quality of 
fit of the function to the experimental data. It may be difficult to 
choose an optimum function from two functions having almost equal 


variance. In this case, the checking of material balance, i.e. the 
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pata 
composition data regenerated from the fitted function, may help one to 
select a function which is consistent with the reaction system. A 
practical case encountered in this work will be presented in the later 
section. As a rule, the material balance is always checked after optimum 
functions for all reacting components have been obtained. 

The superiority of using the smoothing technique of minimizing 
the residual sum of squares using Rosenbrock's Hill-Climbing method 
over the conventional least squares polynomial approximations will be 


presented in chapter 4. 


3.5. Data Smoothing by Minimizing the Residual Sum of Squares and 


Differentiating the Smooth Function 


The procedure for differentiating of experimental data to be 

used in this work is now summarized as shown below: 

STEP 1. A linear or non-linear function, which essentially has 
similar shape to that of the time-composition data is 
fitted to the data. 

STEP 2. The parameter(s) of the fitted function are determined 
by minimizing the residual sum of squares using 
Rosenbrock's Hill-Climbing method (section 3.3.2.) 

STEP 3. The resulting smooth function thus obtained is analytic 


-ally or numerically differentiated. 
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CHAPTER IV 


A NON-LINEAR HYPOTHETICAL: KINETIC MODEL 
4.1. Generating Kinetic Data 
* To examine the accuracy of the: numerical differentiation techni- 
que by minimizing the residual sum of squares using Rosenbrock's Hill- 
Climbing method (section 3.5.), and: the kinetic analysis approach as 
developed in section (2.2.) for complex reaction systems,a non-linear 
hypothetical chemical reaction system, whose rate parameters are known 


a priori, will be used as a reference: system: 


j=2 A 
4 3 
dat 
2A ————> 
L 2 
j=3 Bin 
2 mera 
Figure 14 


Based on this model, time-composition- data may be generated by 
assigning fixed values to the set of reaction rate constants. It is 
assumed that the stoichiometric coefficient of each reacting component 
is equal to the order of the reaction for that component. The object 
of assuming a non-linear model is to: demonstrate the validity of this 
kinetic approach (section 2.2.) i.e. that Peerdae of reaction of each 
reacting component can be obtained directly from the analysis of the 
data without preassigning any value to the order. Generally, imformation 


on the value of the reaction order is not available to the kineticist 
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(especially for heterogeous catalytic reaction systems). The model 
shown in Figure 14 will also be used to illustrate that the approach of 
section (2.2.) treats the time-composition data equally accurate, 
irrespective of the concentration levels, if sufficient data points are 
available. Mathematically, the reaction rate of each reacting component 


for the reaction system of Figure 14 can be expressed by the equation, 


eee eee 
= ib a 
doslack 
geet 52 
de Oy ~ gg + kg), 
aC, 
roe g ra ja? 
dt ak enh SH 
ees OMe 
dt 44.4 


The pre-assigned values for the rate constants and the initial 


conditions are given below (in self-consistent units): 


kay = 1.4, Koo = 0.09 

(4-2) 
ko =80.25, kay = 0.04 
At time t=0, C)=1.0 mole/litre, C,=C,=C,=C,.=0 (4-3) 


The simultaneous non-linear differential equations (4-1) were 


numerically integrated by Fourth-order Runge-Kutta (or Rung-Kutta- 
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Simpson method). With the interval, h=0.002, the accuracy of the 
integration is well assured, as can be seen from the negligible difference 
between values obtained by analytical and numerical integration for the 


concentration, C The material balance column in Table 1 also provides 


he 
a check since at any specific reaction time the equivalent number of 
moles of component, ieee must remain constant. 


The time-composition data thus obtained were simulated according 


to the procedure described in the following section. 


4.2. Simulation of Experimental Data 


First of all, uniformly distributed random number, Xs are 
generated by using IBM 1130 scientific subroutine RANDU and then they 
are inserted into the following equation (4-4) to calculate random 


number, E which are normally distributed around the true values, U. 


(4-4) 
YN/2 


Here, E approaches the normal distribution about the true (or 
mean) value, U, asymptotically as N approaches infinity. With the given 


true, mean value (u), standard deviation (s) is defined by equation (4-5). 


Ss = eU (4-5) 
Then Y=Es+u (4-6) 


The standard deviation s = ev defined in this fashion, corresponds 
somewhat to the constant per cent error experienced in most procedures 
of chemical analysis and the term Es in equation (4-6) is equivalent to 
the normally distributed random error introduced. For e=0.03, as the 


confidence coefficient, Z=25 at a 95% confidence level, Y=it2s, or, 
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ut+0.06u, i.e. a 6% error in the analysis. Similarly, for e=0.05, a 10% 
error in the analysis will be simulated. In this manner, the data can 
be simulated with accuracy as well as the exact knowledge of the true ~ 
values of each rate constant and reaction order. The simulated concent- 
ration data thus obtained were tabulated in Table 2 against their 


corresponding true values. 


es Comparison of Accuracy of Least Squares Polynomial Smoothing 
Technique with the Method in section 3.5. 


Before the method of minimizing the residual sum of squares 
using Rosenbrock's Hill-Climbing method as described in section (3.5.) 
is applied to the analysis of the simulated data for the assumed 
hypothetical model, the accuracy of the method was checked. A comparison 
with conventional linear-least-squares second, third, and fourth order 
polynomial approximations was also provided. For this purpose, component 
As in Figure 14, was simulated by adding 10% nomally distributed random 
error to the true values. These data are also shown in Table 2. 

In the present method, the function, 

k,T 


ke + k,tanh(k,T) 


was fitted to the data, and the computer output is as shown in Table 3. 
The computed values of linear-least-squares second, third, and fourth 
order polynomial approximations are provided in Table 4,5, and 6; 

The smooth data regenerated from the resulted function were listed 
under the title of "present method" in Table 7. The differences between 
this smooth data and the true values are almost negligible. The smooth 
data regenerated from the linear-least-squares second, third, and fourth 
order polynomials were also recorded in the extreme right three columns. 


The plots of the smooth curves obtained by these two methods against the 
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LEAST SQUARES 2 ORDER POLYNOMIALS 


re COEFFICIENTS ‘OF “TRE “PITTED FUNCT TION “ARE 
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TABLE 5 


- 54a - 


LEAST SQUARES 3 ORDER POLYNOMIALS 


THE COEFFICIENTS OF THE FITTED FUNCTION ARE 


AQ ="0001571867 
Al = 0203722293 
A2 ==0200208672 
A3 = 0690011768 
REGENERATED DATA 
X MEASURED Y OBSERVED Y CALCULATED 
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LEAST SQUARES 4 ORDER POLYNOMIALS 


THE CGEPRICIENTS ORS GHE oF ETGEDeEFUNC BION, ARS 


T ERROR = 


AO =-0001006018 
Al = 0002543355 
A2 ==0200066513 
A3 ==0000021039 
AG = 0200001491 
REGENERATED DATA 
X MEASURED Y OBSERVED Y CALCULATED 
1000000 02015470 02014512 
1500000 0024210 0.025959 
22000090 00037490 0.036701 
26500000 02046070 © 0 s0¢6661 
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true and simulated values can be seen in Figure 15,16, and 17. 

The error caused by second order polynomial-fitting in the 
vicinity of the origin are quite significant, although the approximation 
tends to approach to the true values at a distant from the origin. 

Third and fourth order polynomials approximate the data fairly well near 
the origin but deviate considerably as Poacetan time increases. The 
first derivatives obtained from the above four approximations were also 
tabulated in Table 8 to show how the error, which originated from 
smoothing, can be magnified throuh the process of differentiation. 

The fact that ne least-squares-polynomial approximations do not 
smooth the data accurately may be due to the defect that the polynomial 
employed does not satisfy the non-linear properties of the data. Thus, 
the resulting polynomial, obtained by using the criterion of sum-of-the- 
least-squares, may serve as an approximation to the function but not as 
a first derivative, as pointed out by Ralston (aecrien Pee 

As a result, the method (3.5.) is more useful when accurate first 
derivatives are required. This is shown more clearly when the results 


of analyzing the hypothetical model are considered subsequently. 


4.4. Analysis of Kinetic Model using Differentiation Technique and 


Method of Non-linear Estimation 
The simulated concentration data, with 6% normally distributed 
random error at 95% confidence level around the exact values wu in Table 
2, were regarded as the experimental data to be used in this example. 
The data will be differentiated (as per section 3.5.) and fitted by 


non-linear estimation. 
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4.4.1.° Computational Procedure 
STEP: 1 Smoothing Time-Composition Data 


At the beginning of this step, it is: generally useful to plot 
the time-composition data on graph paper, treating time as abscissa 
and composition as ordinate. Visual: examination of such plots will 
give some guidance in selecting simulating functions for each reacting 
component (see Appendices 1 and 2 and references [11,22]). 

The function which has a shape similar to that of the respective 
reacting component is fitted to the data. The parameters present in the 
fitted function are determined by minimizing the difference between 
calculated and true values of the concentration using Rosenbrock's Hill- 
Climbing method. The computer program was written in a format such that 
four different functions could be fitted to the same set of time-compo- 
sition data during a single run of the program. Descriptions of read-in 
statements, input-data and their corresponding format are provided in 
Appendix 3. 

To select the function which: best: fitted the set of time-compo- 
sition data, the variance calculated for each fitting is examined and 
contrasted in Table 9. 

Data at low concentrations are relatively more difficult to 
smooth accurately than data at high: concentrations. Nevertheless, a 
reasonable accurate fit can be achieved when: enough data are available. 
3? Ais and’ A, 
points were employed. The results: of the hill-climbing search for the 


Reacting components, Ay» A serve as good examples; 19 data 
best parameters by Rosenbrock's method are shown in Table 10, 11, 12, 13, 
14, and 15. The initial gwessesused in the iterative hill-climbing sea- 


rch operations, in all cases, were unity. 
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Table 9 


Various smoothing functions for each reacting component 


Reacting No. of Function fitted Variance 
Component £ it 
| 23 
F a babree (Ones, fig, (-042537) 016760210 
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is independent variable, reaction time 
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A3 COMPONENT 
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PARAMETER ESTIMATED BY ROSENBROCK METHOD 
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In Table 9, the results of: fitting: component, ans serve to show 
how one can build a function which: smooths: the data more accurately by 
building upon the knowledge obtained in the preceding trial. In the 
first fit for A, component as shown: in Table 12, the function bate 
(kT) smooths the data resonably accurate except for the first 4 or 5 
Borate: By introducing an exponential term with two parameters to the 


previous function, the accuracy of the fit to the first 4 or 5 points 


improves considerably. 


STEP 2 Material Balance Check 

When all functions which smooth each of the reacting components 
were obtained, the concentration data were: regenerated from the functions 
and the total concentration at: a specific: time: was evaluated and compared 
with the initial value. If poor agreement resulted, then the reacting 
component with the largest variance is approximated again by adding more 
terms to the same function or by introducing: a completely new function. 

As a practical case, example 3: in the following chapter 5, will 
demonstrate how useful this criterion is, especially when the variance 
of the fit of two functions are about’ the same. For example 1, the mat- 
erial balance check was Reedy eee ty tures As can be seen, the 
total number of mole/litre thus obtained from the functions is about 1% 
more or less than that of the initial values. The accuracy here is well 
assured if one recognize the fact: that any error introduced by the smooth 


function, will be magnified four times: in the case of A, component, two 


4 


times in the cases of Ao Ans sgt Me during the process of material 


| 


balance. This is so, because of the non-linear reaction system assumed 


for this example. 
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STEP 3: Differentiation: of the: Function 

The first derivatives of the: smoothed: function can be obtained 
simply by differentiating the fitted function analytically. If the 
function cannot be differentiated analytically, then the five-point 
formula of Whitaker [40] can be used for numerical differentiation. 
The procedure can be summarized as’ follow: generate five equal-spaced 
concentration data points from the fitted function, in such a way that 
the concentration data corresponding to a specific reaction time at 
which the first derivative is desired, is the middle point among the 
five points generated. 

For example 1, the rate balance is shown in Table 17. The non- 
linearity of the reaction system still exists hence the accuracy of the 


total rate balance in Table 17 is:' well: assured. 


STEP. 4 Computation of Reaction. Step 


From the theory developed: in section: (2.2.1.) the rate balance 
for each reacting: component in the: reaction: system of Figure 14, can be 


written as: 


as 11 
SPECTRO Se perme 
tr, = loo (4-7) 
OMe DEVO? alts 
tS = 20 hy 


Because the Left-Hand-Side terms: of equation (4-7) are known from 
STEP 3, and since reaction steps Ta. Loos and Ty, are defined already, 


by substitution, reaction step Tog can be obtained either from the rate 


/ ta xo (OR x 
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balance for A, or 7 components. Component Aas giving r 


- T53 was chosen: because the higher: concentration values available can 


Pha tes BW yee yo) 


eliminate some of the error originating with smoothing. 


STEP 5 Simultaneous estimation of: rate constant and order 
of reaction 
Since the values of the concentrations of each reacting component 
and the rate values of each reaction steps at a specific reaction time, 


have been obtained it remains to find Ks and ny from equation(4-8) 


"i4 = kc . (4-8) 


Because: of the limited knowledge: available on error variance of 
the experimental data, the logarithmic transformation approach suggested 
by Levenspiel [24] is not recommended here. Instead, the Rosenbrock 
Hill-Climbing method is uutilized: to- fit the power function ky jcp 1. 

The smooth values of concentration data, Ci» calculated from the functions 
for each component A, were employed’ in this step. The results for each 
reaction step were presented in the Tables: 18, 19, 20, and 21. 

In brief, the results may be’ summarized with the corresponding 
true values as shown below in Table 22. 

Using the parameters obtained: from these calculations, the reaction 
path is plotted against the simulated data. The agreement is within the 
experimental error even for the low concentration components, Ay» Ay, Ay 
and Ac. The complete plot of all five reacting components is shown on 


Figure 18. An enlarged plot for the low concentration range of component 


Ans Ag» Ay» and As can be seen if Figure 19. 
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TABLE 18 


SIMULTANEOUS ESTIMATING OF RATE CONSTANT AND ORDER OF 
REACTION BY NON@LINEAR FITTING 


ROTATe 


NOe 


EXAMPLE lo» 


REACTION STEP 11 


RATE = RATE CONSTANT*CONC##*ORDER 


RATE CONSTANT N#ORDER 


Oe800000E 
Oel24150E 
Oel30091E 
0e131688E 
0e133526E 
00¢135256E 
0e«l140186E 
02140182E 


00 
01 
Ol 
Ol 
O01 
Ol 
O01 
Ol 


Oe 140000E 
O«l88698E 
0¢192873E 
Oel93884E 
O«el95276E 
00196434E 
DelVois7et 
Os l99735E 


DEVIATION 


=-00155508EF=03 


Oe204741E=04 
00677555E=-04 
Oe 759214E=04 
O0e808388E=04 
02e789538E=04 
Oe 770837E=04 
0e691861E=-04 
0e667274E=04 
0e556409FE=04 
0e544637E=04 
00537056E=04 
Oe432711E=04 
00447630E=04 
003254237E=-04 
00338424E=04 
0e331271E=-04 
00e339839E=04 
Oe270279E=04 


RATE CONSTANT= 
ORDER OF REACTION#=129973 


CALCULATED RATE 


06250414E 00 
00149020E 00 
009873 77E=01 
0e701959E=01 
00524608E=01 
00406889E=01 
06324770E=01 
00265191E=$0l 
00220667E=01 
00186456E=01 
00159644E=01 
00138237E=-01 
00120832E=-01 
0el06547E=01 
00946542E=02 
00846384E=02 
0¢761312E=02 
00688398E=02 
00625702E=02 


16401820 


Ol 
O1 
ot 
Ol 
Ol 
Ol 
Ol 
Ol 


VARTANCE 

00435619E=03 
Oe810760E=05 
Oe314710E=05 
0e229990E=05 
00«136980E=05 
O«e77T5800E=06 
O0e460000E=-08 
Oe460000E=08 


OBS RATE 

00250570E 00 
00149000E 00 
00986 700E=01 
00701200E=01 
00523800E=01 
00e406100E=01 
00324000E=01 
0»264500E=01 
00220000E=01 
00185900E=01 
00159100E=01 
00137700E=01 
00120400F=01 
00106100E=01 
02943000E=02 
00843000F=02 
0«758000E=02 
0 e685000E=02 
0e623000E=02 


CONCe 
0Oe422170E 
00325560E 
00 264930E 
Oe223330E 
06193030E 
Oc l69970E 
06151830E 
Oel37180E 
OsieaLZ0E 
OellS000E 
Oel06400E 
0e99Q0000E=“01 
00e925500E=-01 
Oe 869000E=01 
Oe 819000E=01 
006 774400E=01 
0e734400E=01 
00«698300E=01 
Oe665700E-01 
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TABLE 22 


values obtained 
parameters by this work true values 


kd Pao) 1.4 
My 1.997 2 
Koo 0.094 0.09 
o> 1.017 at 
Ko3 Qe252 On25 
N43 1.000 1 

4d 0.031 0.04 
Diy 0.398 0.5 


4.5. Conclusions 


From the results of example. 1, the: following conclusions may be 
drawn: 

1. the smoothing technique: used’ in: this example gave accurate 
results and its applicability to other’ problems appears justified; 

2. the method of analysis: for multiple-step reaction systems 
successfully simulated example 1; 

3. low concentrations of reacting components can be analyzed 
as accurately as the higher concentration components in the presence of 


enough data. 
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4.6. Testing Alternative Chemical Reaction-Step Modelswwith 
Experimental Rate Data 


In chemical kinetics research, it is always necessary for the 
.dmvestigator to check the adequacy of fit of the experimental rate data 
to alternative: networks of reaction’ steps: which appear to be feasible. 
The following sequential approach was used in this work to evaluate this 
problem. It relies mainly on two criteria: (1) physical realizability; 
and, (2) statistical justification. The former- criterion implies that 
the parameters cannot take eetnesat connacibie with those from nature, 
and the latter criterion rests on relative values of some statistical 
property such as the variance. 


STEP 1 Reaction rates for: each: reaction step, r.,, were determin- 


ij 


ed from the principle of rate balance for: each- assumed chemical reaction 
model. The physical realizability ofeach model was checked by examin- 
ing the values: of reaction rate for each step in the reaction system. 
STEP. 2). Non-linear estimation was’ employed to obtain values for 
at and a The: physical realizability was: again checked by the values 
of parameters thus obtained. 
STEP 3 The variance (or: goodness: of fit) over the entire range 


of data encompassed was obtained: for: each: kinetic model. The variance 


used may be defined as follow: 


m 2 
Lae 
o* = ( i=l i 
bi 


1/2 Z 
TOT : ore 


As an example of this approach, the shape of the time-composition 
simulated data of example 1 suggested the- following alternative networks 


of reaction steps may be feasible’ ones: 
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It is assumed that chemical insights alone do not provide a basis for 


discriminating between these alternative models. 
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In step 1, the calculation of reaction rate values for each 
reaction step, it was observed that the values of Ei = - Si T 53 of 
model (1) are negative in sign. This suggests that the rate of disapp- 
earance of component, a is much faster than the rate of formation of 


A, from Ay: This cannot be true if no amount of Ly is initially present 
since component, Ans appears as an intermediate. Thus, model (1) can be 
eliminated from the subsequent analysis. 

The results from non-linear estimation of Ks, and ne for the 
remaining three models are tabulated in Table 23. 

Models (3) and (4) require negative values for the order of reac- 
tion in some steps and thus are not considered to be physically realiza- 
ble from the kinetics point of view. Model (2) not only shows the low- 
est value of variance for the combined data but also involves reaction 
steps which are all physically realizable. Model (2) is thus the most 
adequate of the four models. 

Since the values of reaction rate at various compositions for each 
reaction step of an assumed model can be estimated, one can generally 
infer some notions about the order of the kinetics to be expected. For 
example, if oT and Ti4 were linearly proportional, first order kinetics 
will probably be suggested by non-linear estimation. On the other hand, 
Lt C, and a were inversely proportional then a negative order of 


kinetics may be anticipated, as for reaction steps, 13 in model (3), and 


both tO and Cay in model (4). 
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CHAPTER V 


APPLICATION OF THE METHOD DEVELOPED FOR ANALYSIS OF 


MULTIPLE-STEP KINETIC DATA 


To demonstrate the applicability of the method developed for 


decomposing multiple-step kinetic models into multiple single-step models 
(section 2.2.), two published reaction systems [7,9], (1) catalytic 
disproportionation of freons, and, (2) catalytic oxidation of naphthalene 


have been re-examined. 


5.1. Example 2, Catalytic Disproportionation of Freons 


Cavaterra, Fattore, and Giordano [7] have studied the catalytic 


disproportionation of CHCL,F, pure chlorofluoromethane, on a highly 


Z 
j fo) 

fluorinated alumina catalyst at a reaction temperature of 200 C. 

Temperature variation along the reactor bed was kept within + 2°! 


Stoichiometrically, this reaction system can be described by the reversible 


reactions, 


2CHCL,F * CHCLF, + CHCL 


2CHCLF,, ee CHF, te CHCL,F 


They assumed that the catalytic disproportionation could be 
represented by a sequence of intermediate steps involving interaction of © 
each chemical species with catalyst; this interaction could be a halogen 
exchange reaction. From the experimental results and the nature of the 


system, they postulated the following pseudo-first-order kinetic. model. 
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Using an abbreviated notation, the chemical model is, 
ki 


z 
ae UL ar Ae ere U 
1 Oko 2 1 


k 


14 
A, + U 5 A 


il Pou (5-2) 


AB pty “23 


Here, A, =CHCL OF, igs A,=CHF.» A, =CHCL., 
and Uo» UL are fluorinating and chlorinating active centers of the 
catalyst, respectively. 

The corresponding kinetic model which follows first-order kinetics 
is described by equations (5-3): 


1 = ! f 2 ' 


waeerany t ee te, t 
de = KU Cy - Ka 9¥yCo - Ky3UQCp 
(5-3) 


= 1 
Cae yk bce 


k14U,C} 


At the steady state, the concentrations of active centers, Us and U,> are 


constant and consequently, equation (5-3) reduces to, 


dc 
oe A Oe ho 


dt i Pw Cy sf ko 209 


14 


dc, 
ae = 101 — Ko9%2 - Ko 3C) 


(S-2) 


te IOHD= ae ee e Dacschecit oa tases em . 
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In brief, the reaction system can be described by the network of 


steps shown below, 


Figure 20 


Applying the method of analyzing first order systems developed by 
Wei and Prater to the reaction system, they found that the results of 
this treatment were consistent with the experimental kinetic data. A 
chemical model involving halogen exchange intermediate reactions between 
reagents and catalyst was used to explain the reaction kinetics within 


experimental error. 


SoZ Object of Present investigation 


It is the aim of this second example to demonstrate again the 
applicability of the method described in section (2.2.2.1),namely that 

(1) the reversible reaction step, in the absence of knowledge of 
the equilibrium constant, can be estimated by non-linear approximation; 

(2) the order of reaction for each reaction step in such a complex 
reaction system may be estimated from the experimental data directly with- 


out any prior restrictive assumptions in the analysis. 
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The experimental data were tabulated as shown in Table 24. The 
reaction times listed in the table were not published but were made 


available through the courtesy of the authors. 


5.3. Application of the Method in Section (2.1.3.) to the Analysis 


The procedure described in section (2.1.3.), when applied to the 
analysis of catalytic disproportionation of freons, may be summarized 


as shown in the following sections. 


Sa oe Smoothing of the Experimental Data 


The smoothing functions fitted to the experimental data for each 


reacting component are listed in Table 25. 


TABLE 25 


Smoothing Function for Each Reacting Component 


Reacting Smoothing Function Variance 
Component of the fit 
CHCL,F Greyre hae leg, Aaa bo O5o5a107 
CHCLF, (1.389T)/[52.48-48.349e 9°97) 41 105e°9-99°M)] 0, 124x107? 
CHF, ~0. O45ARCTAN(T) 9° 3649057097199 T ag, 280 107) 0.133x10* 
CHCL,, 116441149248, 126Tee 9 0.198x10 > 


The data regenerated by the smoothing functions were checked by a 
material balance calculation of the total number of moles present at 
various times. Since the number of moles must remain constant for 


"unimolecular" stoichiometry, the material balances in Table 26 should 
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always approach the value of unity for the assumed. intial number of moles. 
It is obvious that the error is less than 1%. The units used for reacting 


component were mole fraction. 


5.3.2. Reaction Rate 


The reaction rate for each reacting component were obtained by 
differentiating the smoothing functions and were tabulated in Table 27. 
Except the first and last two points, the error in the rate balance is 
less than + 1%. The experimental composition data of CHF component, as 
shown in Table 24, have not only relatively low mole fraction values, but 
also widely scattered data points. Attempts to use less than five 


parameters to determine the smoothing function for CHF, data points were 


3 
made but the smoothing function thus obtained deviated considerably at 
the beginning and end of the reaction period. The use of an empirical 
function with six parameters succeeded in correlating the data points 
reasonably well. This can be seen in Table 26 under the column with 
heading CHF. 
As discussed in chapter 2, when more than five parameters were 
used to characterize the smoothing function, the smoothing tended to 
approach the case of exact polynomials fitting, i.e. although the empirical 
function obtained with a large number of parameters tends to pass through 
most of the data points, only a small amount of smoothing has been done 
by the resulting function. The first value for the reaction rate of 
CHF, at 0.3 second reaction time, was not found to be positive, as it 
should be from kinetic point-of-view, since it only appears as an inter- 


mediate component. This may be attributed to the use of the large number of 
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parameters just described. Because of the erroneous sign for the react- 
ion rate value for CHF. at 0.3 sec reaction time, the rate balance is 
about -0.09 in error, as have shown in Table 27. 

Using the present smoothing: technique, the search for an empirical 
function with less than five parameters to smooth a set of limited and 
widely scattered (+10% error, within 95% confidence level) data points, 
is very difficuly to achieve. This situation improved considerably when 
reasonable numbers of data points were. available, like in example 1 (19 


points). 


5.3.3. Estimation of Reversible Reaction Steps 


As the reacting components A, and Ay, are the two major ones in 


the reaction system, better accuracy in. data smoothing and reaction rate 


for these two may be expected. Because.of this, the reaction steps see 


around the reacting component.A, were chosen to estimate 
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rate constants and order of reactions. for. the forward and backward rever- 


The differential rate. equation for component 


sible steps, ae and loos 


Al» and which describe the kinetic. model will be, 


dC 
1 nj np, --9. . 9022 
Fe yO = Kye, + Ky 0ly 
dc dc n n 
1 4 hwy 11 
i ry ° Sa aE = k GC. = = k C 5-5 
ee) aie eat 22°2 fy >? 


The values of 41 + dC4 calculated from Table 27 and the corre- 
at. at 


sponding smooth values of Cc, and Cc) are listed in Table 20; 
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TABLE 28 


Reaction. Rate i hukor’ and Concentration. of Component Ay» A, 


in Example 2 


Reaction 

tanto. (sec) ecett ah: Wap taeetag sc Aw Ce aces aa) 
0.3 -0.2441 0.0693 0.8575 
0.5 -0.209 0.1061 0.7692 
LO? -0.0812 0.2358 0.4296 
240 -0.0643 0.2526 0.3776 
Pepe -0.0509 0.2653 0.3336 
Zed -0.0374 G2 7.67. 0829 
363 -0.0236 0.2839 O2228 
ae7 -0.0172 052832 0.198 
4.0 -0.0135 0.2801 0.1786 
6.7 -0.0009 O21918 0.0756 
8.0 -0.0032 0.1374 OBO. 


Using. non-linear estimation with Rosenbrock's Hill-Climbing 


method, the values, k 97011648, N, =1.029, k, 70.2964, and Bie + 16% 


were calculated . The remaining two single reaction steps, ThA and 53 


can be estimated by using the following two equations, 
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dC3 dC 
respectively. The values of a aes were listed in Table 27, and the 
corresponding smoothed composition data, C, and Ci» were in Table 26. 


Table 29a gives a comparison of rate constant values obtained by 
this work with those from the relative rate constant matrix, reported 
by Cavaterra et al who employed the methods of kinetic structural 


analysis developed by Wei and Prater 


TABLE 29a 


Comparison of Parameters (relative) obtained by this work with 


Cavaterra's Results 


Cavaterra et al 


This Work 
(Wei-Prater approach) 
Relative Rate order of Absolute Relative Rate order of 
Constant k"’. Reaction n, Rate Constant Constant k'', Reaction n,, - 
aA! sey k akg 14 
ij 
0.2974 at 1.168 
0.2655 0.8927 1.028 
0.1097 0.3689 1.024 
0.0689 0.2316 On 7E2 


The relative rate constants, Ky? obtained from this work and lis- 
ted in Table 28, were obtained by dividing the absolute value of each rate 


constant ke. by that of k 


11° From Table 29a, a meaningfutLocomparisom ‘of ‘ithe 
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relative rate constants reported by Cavaterra with the ones obtained by 
this work cannot be made unless the non-linearity of the order of reaction 
obtained in the present work is taken into consideration. 

A better comparison of the rate constants, taking order of reaction 
into consideration, can be achieved by plotting the reaction paths calculated 
from the kinetic model for both of the cases against the experimental data 
points. To do this, the relative rate constants a obtained by Cavaterra, 
transformed to the absolute values by regarding the rate constant k, ,=0.2655 


14 
for both of the cases. This is justified if we recognise that My of this 
work is approximately first order kinetics. Now, the balance of the rate 
constants can be determined by using relative ratio multiply by the absolute 
value of ke 470-2655. 


TABLE 29b 


Comparison of Parameters (absolute) Obtained by This Work 


with Cavaterra's Results 


Cavaterra et al 


(Wei-Prater approach) This Work 
ij 
Rate Constant Order of Rate Constant Order of 
ke. Reaction n,., ke Reaction n,, 
i af Be i 
at 0.2983 a 0.2974 1.168 
14 Oe2655 a; 0.2655 1.028 
7H 3 21651 . i 0.1097 O24: 


23 OeiZ19 1 0.0689 0.742 
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Using the parameters tabulated in Table 29b and the following 


differential rate equations, to describe: the kinetic model, 


dc n n n 
ae ia 14 29 
a Rie i a 
dc n n n 
2 1 
Ko) a Pecee é AC 
dt (5-7) 
dC; N53 
aha 
dt 
dC Nyy, 
S Rd 
dt 


The reaction paths for both of the cases were calculated. They were 


plotted against the experimental data points in Figure 21. 


5.3.4. Discussions 

In Figure 21, the solid line is the reaction path calculated from 
the kinetic model with rate constants and order of reaction obtained by 
this work. The dashed line represents the composition data calculated. 
for the same kinetic model but using the rate constants listed in Table 
29 under the heading Cavaterra et al, (Wei-Prater approach). 

The calculated reaction paths of both cases for components, CHCLoF 
and CHCL, fit the experimental data equally weil. As for CHF, and CHCL. 
components, the solid lines seem to pass through the characteristic trend 
of most of the data points but fail to predict the last ones, correspon- 


ding to reaction times of 8 secs. This is due to the fact that the smoo- 


thing technique and the non-linear estimation used in this work corresp- 


Pe. 


a _ 
7 : = . 
P| ; 


- 12 + 


i 
P é " ; aes rey 7 
siwollo? oda base d®° sidaT at bataluds? axstemetaq ofa g ra _ 


ad 


: oe °. . A eT Ae 
~labom sitentiA edz ditsesb o3 ,anolssups ois1 Istsasze22th 


e oe ec } si a ’ «Ub » 
nae | - gIoe4 7; 7 3 a a 
(T=2) | ais Paid P onvdin 


28880 93 Yo diod +03 eaddeq aoljosex 


it 
: 
t 
T 
4 
WZ 
? 
its 
<4 
= 
FO 
nr 


Sisw ysqdT .bs 


c + t 


* 
IS sigtt of e3 sata siab {eanemtregxs ort Jentsgs besz0! iq 


= * o 


- 7 ’ 


mott bs3sivolss dteq noltsosex sd eb enti biloa sz AS stugtl ot © 
an 


“3 - an * 
rd bante3do noljosst to srebto bab a3 a8: B09 o26t. iste febom otteah f 


bsisluolas sisb nolialecumos end etn aesiqer past endead od? wb 


“ 


mt bsjatl einssenos: 448% $63 ‘soba aud Lsbom ‘thoatia onan od 
. (d5s01qqs t93819-IoW) ,is Je, exzotevad gnibeod. os ; sbau eS 


me 
GoJDHD .s3nsnoqmon yor as s> riod 26 eridaq nohvanen besatuotaa oa 


g40HO bas THO rot eA “Ateu vilsups sieb Jadnsetaogne ony Pre 


‘te an 

brova aiteiiss opis on ds®addy eesq of wa9e pubennin | 
a ' \ 

oqe ItOD ,e800 desl siz Jotbetq ot tre a be iosaen Te 

ae ee * t 


1 y ¢ = 


va 
9. 308 » 92. 02. sub ate i 


a | + 
& . 9 7 
f ai beey cobfen , > 


S SOON SILANIM 30 NOSTAVGNOD $12 ean3tg 


@2ie)- SN NOV 


OR" i 


HIvOdddY HAlVMe - Da 
ONVONOTS SVoealVAV = SNI1 CGHSvo 
SOOM STHL = GNI OF 


eset = CSOBDI EYE-LET SG MRIS Bae A 


edhe & “iSCWiVd oe Br Ctclekal helene s 
S 1H 6 on SUNT 
“wos ohHS. % TWINING 


NOLLNOTIM EAST ILA WLVI = WaALSAS 


TION 


NOLTLIVaS 


I =e a : S 2.” 
5 THs: ei 2) an oS TA. 
va 1 eer rh J OW" 
| SP tl = wld! 4 


= 
~ 
“ 
— 
D — 
, , 
- - 
—— — fa a} 
rect s 


= 99. — 


ond somewhat to equal weight for all data points. Here, although it is 
difficult to conclude which reaction path predicts the experimental data 
more satisfactorily, the following conclusions may be drawn for the 
method used in this work. 

(1) Order of reaction for each reaction step can be estimated from 
the experimental data directly without any restrictive assumptions such 
as all reaction steps follow first order kinetics that Cavaterra and his 
co-workers have made in their analysis. 

(2) The reversible reaction step can be approximated by non-linear 
estimation as described in section (2.2.2.1.) without the equilibrium 


constant. 
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4 Example 3,.Catalytic Oxidation. of. Naphthalene 


As desribed earlier in section (2.2.2.), the apparent reaction steps 
in this reaction system formed a closed loop. This example demonstrates the 
applicability of the method described in section (2.2.2.). In view of the 
importance of the catalytic oxidation of naphthalene reaction in the manu- 
facture of phthalic anhydride, further study is.of interest. The data of 
D'Alessandro and Farkas [ 9] were employed.in this example. It is also 
worth noting that this is a multiple-step heterogeneous catalytic reaction 
system. Irrespective of diffusional. problems. which may exist the model- 
ling procedure is analogous to that employed. for homogeneous reaction 
systems. 

D'Alessandro and Farkas studied the vapor. phase catalytic. oxidation 
of naphthalene over a catalyst of. vanadium pentoxide. The reaction was 
carried out in the presence of a large excess. of. oxygen under conditions 
such that a negligible change in volume. was encountered over a temperature 
range from 340°C to 475°C. 

In their mechanistic analysis, they first. grouped the highly oxidi- 
zed and relatively low conversion products like maleic anhydride, canhers 
monoxide, and carbon dioxide. They. proposed. the. reaction scheme of Figure 
22 with the restrictive assumption that all reaction steps follow first 
order kinetics. Using the D'Alessandro's abbreviation for the reacting 


components the chemical model is as shown in Figure 22. 
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Figure 228 


The differential rate equations for this. system will be, 


dC 
per Ne ey 
gee 7 Oa? Rips I IG. 


de *11%Np ~ (K24.+ 26) Cnq 
(o=8) 


dt Ko4Cq + k12oNp + k35Cpa 


at £13°Np + K26nq + *35°pa 


Because the solution to the above system of differential ‘equations 
cannot be obtained explicitly, they suggested. obtaining the. relative rate 
constants of PA, Nq, CO, first by extrapolation of the product distribution 
to zero conversion. Since these are assumed to be. first order plots, the 
relative values of the initial rates are proportional to the intercepts 
thus obtained. If it is assumed that the rate constants are also proport- 
tonal to the initial rates (which is a mathematical consequence of the 
assumed linear kinetics), then one obtains the ratios, 
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The remaining relative values of the rate. constants, Koy ko6s and 
k35 are adjusted to the proposed model by trial-and-error. The data used 
in their analysis was measured at the. reaction. temperature of 410 C. 
Two different sets of values for the unknown rate constants were reported, 
and a comparison of the composition distribution curves plotted from these 
models, with that of the experimental ones was made. They found neither 
model satisfied the experimental data. As a remedy, they argue that since 
the concentration of carbon monoxide plue carbon dioxide does not increase 
with conversion, whereas the concentration of maleic anhydride does increase 


with conversion, a better reaction scheme will be given by 


PA “34 MA 
——————> 
k11 
k93 
ky 
° 
Np 
5 
CO + CO, 
Figure 23 


In this thesis, the latter. chemical. model. is. called: model 2. 


5.4.1. Smoothing and Differentiation of. the Experimental Data 


Because of the limited amount. of data available at reaction tempera- 
fe) 
tures other than 410 C, and the scattering. apparent in the limited data, in 


this work, only the data at 410°C were examined. The isothermal experimen- 


tal data exployed are shown in Table 30. 


An illustration of how the criterion of a material balance check can 
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help one to choose the "smoothing fitting" from two different fitted 
functions, whose variances are about the same, will be presented. Table 
31 shows a material balance which was checked by regenerating the 
composition data from the two fitted functions. The phthalic anhydride 
reacting component was used for the material balance and the smoothing 


functions employed were 


(a) 2.0302TANH(T) 2°73 and 


(b) 0.7563TANH(5.67T) 


From Table 32, it is obvious that the better material balance for 
the whole reaction system was obtained fitting the second of the two 
functions. Though their error variances are approximately equal, the 
former being 0.123 x WO; sbandethe dacter. 02139A x eee the second function 
was chosen because of its better material balance. 

In calculating the compositions of Table 32, the functions used 
to fix the data for the balance of the reacting components are shown in 


Table. 33, 
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Smoothing Functions for Each Reacting Component (Example 3) 
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The reaction time, T, was used as the independent variable in the 
fitted function. Because of the wide scattering of the experimental data 
points (especially for the low concentration values for reacting components 
like maleic anhydride and carbon monoxide, and the limited number of time- 
composition data points available, the material balance is not very 
satisfactory for the first point in Table 31. The wide scattering of the 
data may be due to a 7°¢ range in catalyst temperature during the 
experiments. 

Table 34 shows the first derivatives obtained by differentiating 
the fitted functions corresponding to each of the six reacting components. 
The sum of reaction rates for all components at a common instantaneous 
time must equal zero. This "rate balance" shown in the extreme right 
column is not very satisfactory, especially the first point. The same 
reasons for the inconsistent material balance just mentioned are applicable 
and the error originated from the fitted functions is magnified considerably . 
by taking the first derivatives. 

At this stage of the calculations, all of the reaction rate values 
for the reacting components and the corresponding smoothed composition 
data are known. In addition to the kinetic model 2, model 1, proposed by 
Ioffe et al [18] during their study of the same reaction system as described 


earlier in section (1.1.2.) is shown below 
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5.4.2. Selection of the Best-Fitting Chemical Model 


Chemical model 1 


Se 
k25 
Nq ————» CO + CO.. 


Chemical model 2 


Applying the procedure for the loop. reaction system described in 


section (2.2.2.), reaction steps rj, and. ro, were.estimated. by non-linear 


estimation. 


From the principle of rate balance, both models give 
Toa = Ty1 a To3 - Tya> 


~ nd neo . 


The reasons for choosing reaction steps T11 and 93 around PA 


reacting component are as follows: 


PA is a major component in the reaction system and,.as a result, 


(1) 
the first derivative can be obtained more accurately than for the 
rest of the components; 

(2) the reaction steps, rj, and Yo3> are equivalent in both the 


models, i.e. the values of kiq N11? ko3 and Ny3 obtained from 
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non-linear estimation apply to both models. 


From Table 34 the left-hand-side of-equation (5-9), r 7 can 


PA *' TMA 
be calculated. The rate data thus obtained and the corresponding concen- 


trations of Cypand CN are listed in Table 35. 


Table 35 


Values of reaction rate (rp, + Ta) and concentration Cy? Cyg 


Reaction time (sec.) 


In equation (5-9), k ko3 and N53 are unknown parameters 


ee? 
which may be estimated by fitting the function, ky (Cy) “11 a ky3 (Cy) "23 
to the rate data and by using Rosenbrock's Hill-climbing Method to minimize 
the residual sum of the squares between the observed. rate values and that 
of the calculated ones. The parameters. were. found to be optional when 


Ue ante I= SSE ig 0.7154, ko3 = 1.2298, and no3 = 0.8935, with a variance 


of 0.0354. 


For Kinetic model 1, since the rates for the reaction step, TH 
n 
) 


can be calculated from, rj, = ki (Cyp 12 


Aig 3k se 
seLnUuse f “Np t4, can also 
be calculated. The rate values for reaction steps, l49° To5» and fz, are 


recorded in Table 36. 
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Table 36 


Values of rate for reaction step, T19>9 93> and T3, 


Reaction time (sec.) 


13112 


"30.065 0.9436 0.487 
0.086 0.6498 0.4158 
0.13 0.237 0.2971 
Onl7e -0.2695 0.2163 

-0.1818 


In the same manner, the reaction step, Tyo» in Kinetic model 2 can 


be calculated from PY>, = “Typ =iFy) 7 T15- 


Table 37 


Values of rate for reaction steps, Tio» T34> and ta 


etm: Sie aide Bd 4 oe a_i Le Cea 


~ FCOo + CO 


Reaction step 3, in both of the models is a reaction between the interme- 


diate PA and the final product MA, and so the estimated k 4 and n_,, values 


3 34 


are good for both reaction schemes. 


Routine application of the procedure of non-linear estimation for 


each reaction step yields the results tabulated in Table 37. 


The kinetic models for the two chemical schemes are given by equa- 
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tions (5-10) and (5-11). For model 1, 
dc 
Np “big n 
ae 7 wk11Snp  - K12Cyp 22 
dc 
PA n B 
f (5-10) 
ong =4 my P53 N95 
= Cc See = koec 
dc 
dt 34°PA 
“Sch, + co 1), 25 
dt 25°Nq 
and for model 2 
dc. 3 
Dig = ss ie 1 
dt Kron Kol, * = K15¢np . 
dC, 4 : : 
dt rate 23°n Sey 
ei (5-11) 
Na — n . 
= 12 
ie 12H ¥93°N 
dc 
MA n 
k 34 
dt 34°p a 
dc 
CO, + CO i 
dt Kky5Cy, 19 


Integration of these kinetic models was possible by using the parameter 


values listed in Table 38 and a fourth order Runge-Kutta method. 


The reaction paths calculated from the Kinetic models were 
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compared with the experimental data. They are as shown in figures 24 


and 25. 


In figure 24, the calculated reaction paths are in good agreement 
with the experimental data, except for that of maleic anhydride which is 
slightly lower than the experimental one...On. the whole, the. accuracy is 
within the scattering of the experimental data. In figure 25, the calcula- 
ted path of PA is in good agreement but the rest of the reacting components 
deviate significantly from the experimental data. From the results of the 
calculated reactdon paths and the variances for the combined data, it 
appears that Kinetic model 1 fits the given experimental data more satis- 


factorily than model 2 proposed by D'Alessandro and Farkas. 


Using the data of Ioffe and Sherman [18], Peterson: [29] has investi- 
gated the same vapor phase oxidation of Naphthalene by non-linear estimation 
for the data between 260-400'C reaction temperature. The- approach he used, 
of the type in category (2) in chapter 1, involves the numerical integrations 
of the differential equations and subsequently applying non-linear least 
squares to fit the data. From his analysis, he found that both Kinetic. 


models 1 and 2, fitted the data equally well. 


In Table 38, the values shown for orders of reaction, n3, "= 2.485 
and Mio ® 2.8279, are physically meaningless. This is perhaps: to be ex- 
pected since the kinetic models are not comparable in form to mechanis- 
tic models such as Langmuir-Hinshelwood models. In addition, the use of 
six data points considerably restricts the validity of inferences which 
may be drawn. iduartaleee the use of combined single-step power-law 


rate equations in modelling of complex reaction systems introduces an 
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element of compatibility between chemical and kinetic models. This as- 
pects is missing when a complete model is assumed, and not constructed 


from the individual steps. 
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CHAPTER VI 


DISCUSSION 


Ou, Summary of the Method 


The method proposed in this thesis for the interpretation of experi- 


mental isothermal kinetic data for multiple-step chemical reaction systems 


can be summarized as follows: 


Sree Ls 


STEP II. Develo ment of Alternative Chemical Reaction Ste .Models 


STEP III.Decomposition of Overall Rate. Data Into Single Step: Rate Data 


Smoothing and Differentiation of the Experimental Data 


A linear or non-linear function of similar shape is fitted . The 
parameters of the fitted function are estimated by minimizing 

the residual sum of the squares of differences between calcula- 
ted and observed compositions using Rosenbrock's Hill-Climbing 
Method. Differentiation of the resulting functions provides local 


values of the reaction rates for each reacting component. 


If the overall chemical steps for the reaction system being 
calculated are unknown, networks of apparently feasible reaction 
steps must be devised. Such alternative chemical models are 
based upon the sequential and parallel chemical steps deduced 
from the eet deed ia ascii plots and/or additional chemical 


insights. 


Decompose the overall rate data for each reacting species into 
rate data for individual reaction steps using the rate balance me- 
thod. A non-linear approximation for the estimation of two react- 


ion steps is required when the coefficient matrix of the reaction 
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system is singular. 


STEP IV. Evaluation of Kinetic Parameters for each Reacting Step in the 


Kinetic Model 
The parameters, Big and ck’, of the power function rate equation 


for each reaction step are estimated by minimizing the residual 


sum of squares as for step I. 


STREP OV. Discriminating between Chemical Models 


Alternative chemical reaction-step models are checked by a 
sequential approach which relies on two criteria 
(1) physical realizability, and 


(2) statistical justification 


6.2. Comments on Applicability of this Modelling Method 


The procecure of Step I is obviously restricted to suitable time- 
composition data such as that measured with batch or integral reactors. 
However, at this stage of development, the application of the method to 
homogeneous system or to catalytic heterogeneous systems is identical in 
procedure. If differential rate data are available directly from the 
experiments, they may require smoothing but not differentiation. The 
method is not recommended for use when limited experimental data are 
available because of the resultant uncertainties inherent in curve-fitting 
methods. Five or six data points were considered to be minimal for 
modelling depending upon the complexity of the chemical models involved. 
The relevance of the. kinetic model to the true chemical model is questio- 
nable when the time-concentration profiles for various components are 
clearly ambiguous, i.e. when the scatter in the data necessitates 


approximating of unknown functional forms. 
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When the requirements of. sufficient. kinetic data and clearly de- 
fined concentration-time trajectories are met, this approach to kinetic 
modelling should prove to be superior to other methods of analyzing non-. 
linear models which model the complete. kinetics. simultaneously. This 
assertion is made on the basis that the method discriminates between the 
non-linearities of individual steps in. the chemical model: before cons- 
tructing the overall kinetic model. With prior integral methods of 
analysis, only trial-and-error can. be used in constructing suitable models, 
With complex systems, mathematical approximations rather than. chemically 


relevant kinetics models result. 


6.3. Conclusions 

1. The first derivatives obtained by differentiating: the smoothing 
functions. used in this lee a aceivere venules and their. applicability 
to the interpretation of rate data for.a multiple-step reaction system 


appears justified. 


2. The order of reaction and rate constant for each reaction step 
can be evaluated from the isothermal experimental data without the con- 
ventional trial-and-error modelling procedure often used in. other non- 


linear estimation methods. 


3. The use of both integral and. smoothed differential: rate data 
provides an improved basis for discriminating between alternative chemical 


models prior. to development of the kinetic model. 


4. Situations involving reversible or closed (loop) reaction steps in 


a multiple-step chemical reaction system pose computational. difficulties 
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in this method. This has been successfully. resolved by use: of. non-linear 


approximations for these steps before. decomposition of the.network of steps. 


5. The iterative computational procedure involved is-very simple. The 
use of a direct search optimization procedure in evaluating parameters for 
each step should pose no problems since only two to five- parameters are 


generally required per step. 


6.4. Recommendations 

1. The analysis of error variance of experimental data a priori to the 
application of the method is desirable. The correlation between the error 
variance of the experimental data with the accuracy of the: parameters estir 


mated by the method should be examined in more detail. 


2. Sets of time-composition data with different initial values. are re- 
commended to obtain more accurate and consistent values of: the- parameters. 
This may be particularly effective. in determining the value- of-the order of 


reaction for each component. 


3. Use of spline function to smooth a limited experimental data 


available is worth considering. 
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NOMENCLATURE 


denote a column matrix (ay, Any teeeeeseeceeeceny ay) 

parameter of the approximation to the rate equations, dimensionless 
the vector of the unknown coefficients, Ba defined in equation (3.8) 
concentration of reacting component A at time = t 

initial concentration of A when time = 0 

concentration of qth component 

error constant factor 

the error which is assumed to be randomly distributed, i.e. with 
mean uw = 0, variance 52 

normally distributed error around the true value, u 

concentration dependent functional relationship terms of the rate 
equation, defined in equation. (1-9) 

molal feed rate, lb-moles. feed/hr. in. equation (3-14) 

interval for equal spacing 

reaction rate constant 


th reaction 


dimensionless rate constant for the i 
reaction rate constant: in- equation. (5-3) 

rate constant in bareteien (1-11) 

the degree of polynomial, defined in equation (3-4) 

the number of reacting componets in the reaction system in equation 
(4-9) 

coefficient matrix of the reaction step 


th component in ach reaction 


order of reaction 1nepe ab 
number of values of Xs to be used 


number of data points in equation (3-17) 


number of parameter in the fitted function in equation (3-17) 
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P(X; ) some as yet unspecified function of. independent variable x and 
is of degree j 
Q the number of time intervals between the time at which the 


concentrations were measured 


; eth 
at _ reaction rate for i+ component 
> 
= Eas 
j=1 1j 
: Aaa weet : 
Ti reaction rate for i component takes part in j reaction 


= f(t, Cy» Cos Cas « SUS. PF CODE=HN, BPLyHomien) pa( poseLpe2, (Ir45-/n) 


ry column matrix of reaction rate for oP component (ry Tos re a) 
Ys. column matrix of reaction step r,, 
ij ij 
s standard deviation in equation(4.5) 
t reaction time 
a independent variable, reaction. time 

= 
ig the time required to reach 1 for.an. arbitrary initial value of 8 

C 
tl the time required to reach a given.value of y (= ) when 
Ao 

By = ] 
u thesindex:*invequation G11) (a= 1, 2, sie og, Q) 
Ke fluorinating active centers of the catalyst 
v1 chlorinating active centers of the. catalyst 
W mass of catalyst, lb 
x conversion 

S Oe ok 
x a uniformly distributed random number 
X denote a vector of column. matrix (X15 Xo» sores ets X,) 


4 represents the resultant simulated data 
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dimensionless 

a Cy/Cyy 

Cy(ty) - C,(tQ) 

experimental concentration values of. reacting component 

the 4 polynomial coefficient, fefined. in equation (3-4) 

the coefficients of the power series 

— dimensionless, subscript r.denotes a reference ca ge ees 
the coefficients of the orthogonal polynomial in equation (3-11) 
Variance for the complete an of data 

the variance of the fit for as component 

mean value 


determinant multiresponse error_ criterion 


integration constant in equation (1-26) 
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SOME COMMON FUNCTIONS. THAT. TAKE. THE FORMS 


OF THE REACTANT DISTRIBUTION. KINETIC DATA 


eMAOT BHT 


ATAC O14 


radIA MOITUATATEIO TMATOARA aut YO 


f-L 


fe 


fi 


, 


oe hed? 


‘ft MTOMaTEA - 


~ 


aye 6°) AST OTA * ‘ a ¢, 7 ? _ ie toy! : 35 
\T TART 2yOTTOMUT MOMMOO-amoe® "" - n 


at 


£ 


1-2 


Beep OTZSUTH uvoTAnqTaqstTp Jueqzoerer |yq jo 


QO a, @ a) QO 
< c e = & 
S nm a _ as 


SUIOJ 9Y4} aye EY suoTRoUNG 
X 

O i O O 

ie gt a ae 


Ante Oona Ort Aa te nae UA 


(AeS*O-IgGxa=A 


PAC bt bd 7 be Ate 


uowulodD 


iO 


awosS T[ xtpueddy 


CMMI O 
SAAS 


OOOT + 


0002 * O° 


GEE" -O 
(A 
BEES °C 
GEES<C 
BEE/°C 
66a8+0 


Estate; ¢O 


‘2 =F 


“ees ot3* si=¥ F 


i%, 
/ : | 


(Wal <~ 1 SK a=¥ al a). 
- } | 


2-1 


APPENDIX 2 


SOME COMMON FUNCTIONS. THAT. TAKE. THE FORMS 


OF THE PRODUCT DISTRIBUTION. KINETIC DATA 
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DESCRIPTION OF INPUT. AND. READ. IN. FORMAT FOR 


ROSENBROCK!S METHOD 


cs 
* bo 
~ Bae 
- / 


AO% TAMAOT M1 GASH GMA TUIMI YO WOTTaIADesa 
H a 7 I 
a. ey: : ; i ~ee 


cOHTaM 2'wposéumzon = .. , 


< 
€ 


3-2 


Appendix 3 


Rosenbrock"s Method 


Description of Input Data and the read in Format for Rosenbrock's Hill- 
Climbing Method. 


M number of variables 


N 


number of data points 
NFIT = number of fits 
p = number of parameters 


KPN = O if normalization desired 


Imi not 


M1 


maximum number of rotation (=50*P) 

M2 = maximum number of steps between rotation (=75%*P) 
FORMAT © (615) 

XNAM(I,J) = title for the independent variable (Time) 

DES(J) = description of the input function 
“FORMAT (20A3) 

Y(1) = dependent variable (concentration for the present case) 
FORMAT (F13.7) 

X(I, 1) = independent variable (Reaction time) 

ss = a constant 

PAR1(2) = Parameter (I = 1,p) 


FORMAT (6F12.7) 
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APPENDIX 4 = (1 
c ROSEN METHOD 
ROSEN METHOD 


HHH HH HHH HIE HSE HE HE FEE HEHE HHH HEHE HHH HH HH 


# * 
% ROSENBROCK'S HILL@CLIMBING METHOD * 
# 


FE Fe HH FE HE HEE FE HT FE FE FE HE HE HE SE FE HE HE Fe HE HEHE FE TE HE HE HE HE HE FE HE FE HE He HE HE HE 


INTEGER P 

DIMENSION ELAM(9) sDES(81) »XNAM(999) 9KS(9) 9O0V(999) 9X(40 
#91) 9Y(40) oPA ; 
1R1(9) sPAR3(9) sKSF(9) KF (9) 


THIS IS THE MAINLINE PROGRAM FOR ROSENBROCK'S METHOD 


INPUT DATA READIN 


MXMN=O 

READ (591) MoN 
READ(592) ( (XNAM( 19J) 9J=195) sIT=1l 9M) 

RED $592) NEFITS 

DO 3 T=loeN 

READ oS) ACL 8 

CONTINUE 

DO 4 T=l9N 

READ (595) Y(1I) 

CONTINUE 

FORMAT (Fl3e7) 

DO 6 NFT=leNFITS 

READ (591) PsKPNoM1l M2 

FORMAT (615) 

READ(592) (DES(J) sJ=1924) 

FORMAT (24A3) 

KP=P+1] 

READ. (S97) Sost PARITY) sl=lsP) 

FORMAT (9F8e4) 

WRITE (698) 
FORMAT('1's///////948X9'AB COMPONENT! 9/948X9'e= eeeen= 
tmanl) 

WRITE (699) (DES(JU) 9y=1924) 

FORMAT ('0'935Xs24A39/) 

WRITE (6910) 

FORMAT (' '932X9'PARAMETER ESTIMATES BY ROSENBROCK 
* METHOD! 9/) 


GENERATE INITIAL SET OF ORTHOGONAL VECTORS & NULL 
SCALING FACTORS 
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C ROSEN = METHOD eee (CONT!D) | 


DO 11 I=loP 
KSF(I)=0 
PAR3(1I)=PAR1(1) 
DO 12 T=l9P 

DO =P3- VElsP 
OV(I9U)=00e0 
OV(I9T)=1¢e0 


ADJUSTMENT OF PARAMETERS 


CALL VARTA (PAR1 sKSFePsVAR2  oNoNFT oXeY) 
VARI=EVAR2* (=1le) **MXMN 

VAR3=VARI1 

COUNT=06 

DO 14 KKK=1l0M™1 

DO 15 JT=loeP 

ELAM(1)=SS*(100) ¥*KSF(I) 

KS(1I)=0 

KF(1T)=0 

DO 16 KLP=19M2 

DO 17 ITIT=1l9P 

DO 18 I=loeP 
PARI(I)=S=PARI(1)+OV( II 9I)*ELAM(IT) 
CALL VARIA (PAR1 6KSFoP9VAR2Z oNoNFT 9X9) 
VAR2Z=VAR2*(=le) #¥#MXMN 
COUNT=COUNT+1le 

IF (VARI=VAR2) 19920920 

VAR1=VAR2 

KS(II)=1 

ELAM(II)=ELAM(II)*320 

GO s10e 24 

DO 22 T=leP 
PARI(I)=PARI(I)@ELAM(II)*OVIITI 91) 
CONTINUE 

ELAM(IT)=ELAM(I1I)/(=200) 
KFCII)=KS(IT) 

KS(II)=0 

DO 23 T=1leP 

CT=ELAM( IT) /PARI(1) 

CTEST=ABS(CT) 

IF(CTEST=1leE"08) 24924925 

IF (KF(I)) 24917924 

CONTINUE 

GO-TO 26 

CONTINUE 

CONTINUE 

CONTINUE 


ROTATION OF AXES 
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APPENDIX 4 = 3 
C ROSEN = METHOD eee (CONT'D) 
CALL ROTAT(P»sPAR1 sPAR3 s0V) 


DATA OUTPUT 


DO 27 IT=loP 
PARICI)=PARI(I)/((100) *#*KSF(I)) 
2? CONTINUE 
WRITE (6928) KKK o(PAR1(I)sI1=19P) sVAR1 »COUNT 
28 FORMAT (* '916X9I301X910E1205) 
DO 29 JT=loP 
PARI(T)=PARL(I)*( (100) *#KSF(I)) 
29 PAR2(1)=PAR1(1) 
IF (ABS(VARI=VAR3)-1leE=12) 30930931 
31 VAR3=VARI1 
IF (KPN) 149309014 
14 CONTINUE 
SO LF URPN=1) 3is32 93:1 
31 KPN=1 
CALL TORM (P»9PAR19KSF) 
DO 33 T=leP 
33 PAR3(1)=PARI{(1) 
GO TO 34 
A2° D035 121 9P 
35 PARL(I)=PARI(1I)/( (100) #*#KSF(I)) 
WRITE (6936) 
36 FORMAT ('0'9/s923Xe'DEVIATION! 96X9'CAL FUNC VALUE 
* OBS FUNC VALUE 
1's6Xs "REACTION TIME*) 
DO 37 JFloeN 
FUNC=06 


FOR DATA CORRELATIONs INSERT FUNCTION AFTER THIS CARDe 


GO TO (200920192029203) »sNFT 
200 FUNC=PARI (1) *EXP(PAR1(2)*X (491) )+PAR1(3)*TANH(PAR1(4) 
¥#X(J91)) 
GO TO 210 
201 CONTINUE 
GO 70 210 
202 “CONTINUE 
GO FO. 210 
203 CONTINUE 
210 CONTINUE 
DEV=Y(J)=FUNC 
WRITE (6936) DEVsFUNCSY(40) 9 (X(U9T) 9 T=1l9M) 
326 FORMAT (! 'ol9OXsElLS eB 92Xstloebe3XsE—lceb95XvEl15 06) 
37 CONTINUE 
6 CONTINUE 
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APPENDIX 4 = 
C ROSEN METHOD eee (CONT'D) 


WRITE (6938) PAR1(1) »PAR1 (2) »sPAR1 (3) sPAR1 (4) sVAR1 
38 FORMAT ('O'o//o27Xo*KI=" sFBeS92Xa'K2=!' oF Be592X9'KI=! 
*9FBe592Xe'KGes! 
LesFGe593Xs'!VARIANCE=! 9F1l008) 
S FOR 
END 


4 


o 


-Q XIOWBSSA van 


; “a x 7 
) (O'TAOD) ees GOHTSM waeOR 3 — 


LAAVs (0) DAGe (EWA (SP CRAG SLAG 1BEeS a 
Teen  eXSees ah 4 AMAR sé 
itemena 
(340 ABDHATBAY MEUE SO AEL 
a0Te 


OW3 


i 


(TOY te 


2C0 


201 


202 


203 


fo 


Ww 


ABPENOIX a ie 48 
SUBROUTINE VARIA | 
SUBROUTINE VARIA (PAR1LeKSF sPsVAR2 sNoNFT 9XeY) 


INTEGER P 
DIMENSION PAR1(9) 9KSF(9) 9X(4091) s¥(40) 


THIS SUBROUTINE CALCULATES THE VALUE OF THE FUNCTION 


10 BE 


CPTIMIZEDs ALSO THE VARIANCE IN THE CASE OF DATA 
CORRELATION 


CO 5- l=leP 
PARL(I)=PARICI)/ (0100) #*KSFI(I)) 
DEV2=0e0 

DO 2 J=leN 


FOR DATA CORRELATIONs INSERT FUNCTION AFTER THIS CARDe 


GO TO (200920192029203) »NFT 
FUNC=PARI (1) *¥EXP(PARL(2)*#X (491) )+PAR1(3)*TANH(PARI(4) 


##X (Jel) ) 


G0. TO, 2b0 
CONTINUE 

GO, TO 72.0 
CONTINUE 

GO. 7D, 210 
CONTINUE 

oO TO 1210 
CONTINUE 
DEV2=(Y (J) =FUNC) ##2+DEV2 
CONTINUE 
VAR2Z=DEV2/(N=4) 


DO 3 JT=l9P 
PARL( I) =PARL(I)*( (100) **KSFIIT)) | 
RETURN 


END 
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SUBROUTINE ROTAT 


SUBROUTINE ROTAT (P»PAR1 »PAR390V) 
ENTEGER: P 
DIMENSION A(999) 9B(999) sPAR1(9) sPAR3(9) s0V(999) 


THIS SUBROUTINE ROTATES THE AXES USING THE GRAM 
SCHMIOT 


ORTHOGONALIZATION PROCESS 


DO 1 T=loepP 

DO 1 J=loeP 
A(lsJ)=00e0 
BlIsJ)=0e0 

DO 2 T=loeP 

DO 2 J=IeP 
A(IsJ)=PAR1(J)=PAR3 (UU) 
Bi IsJ)H=A( I 9J) 
SO=0¢e0 

DO 3 J=loP 
SO=SOtFA (lou) #*2 
S=SO* #005 

DO 4 J=loeP 

OV( lsd) =A(19U)/S 
DO 5 IT=29P 
MI=Iel1 

DO 6 K=l»M!I 
DOG=0e0 

DO 7 J=loP 
DOG=DOG+OVI(K9U) FAK oJ) 
DO 8 J=EloP 

Bl IsJ)=B(I194U)=DOG#OVIK 9) 
CONTINUE 

SO=020 

DO 9 J=loeP 
SO=SO0+B( 194) ¥¥*2 
S=SO**0 05 

DO 10 J=lo»sP 
OVOPsJIT=BC hs JI7S 
CONTINUE 

RETURN 

END 
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SUBROUTINE TORM 


SUBROUTINE TORM (PsPAR1»sKSF) 
INTEGER -'P 
DIMENSION PAR1(9) 9KSF(9) 


THIS SUBROUTINE MAPS THE PARAMETERS BETWEEN Ool & 160 


DoT? fails? 
KSF(1T)=0 

DO 2 [=19P 
PP1L=ALOG(PAR1(1))/2e303 
IF(PP1) 39494 
KSF(I)=PP1l=1 
GO.ARO SZ 
KSF(I)=PP1 
CONTINUE 
RETURN 

END 
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Cc RUNGE = KUTTA 
G RUNGE = KUTTA 
C HEH HEE FE HEH HEH FE AE FE HE HE EEE HEHE HE HE HE HE HE HE IE HE HE SE SE EH Ee HH HHH 
G + 
C * . NUMERICAL INTEGRATION BY FOURTH ORDER ¥* 
} * RUNGE=KUTTA'S (OR RUNGE@KUTTA@SIMPSON) * 
C * METHOC * 
c * + 
G He FE HE HEHE SESE EE ESE ESE EEK EH HEHEHE FE SEE TE HHH HHH HH HHH HHH 
C 
fe 
DIMENSION B(5)9C(5) 9D(5) 9E(5) 9CA1(5) 9CA2(5) 9CA3(5) 
¥9CA4(5) 9CA5(5) 9X 
1CA1(5001) »XCA2(5001) 9XCA3B (5001) 9XCA4(5001) 9XCA6 (5001) 
%#97T(5001) 9wl5) | 
19XCA5 (5001) 9XM(5001) 
DX=0el 
M=5 
MRSS 
READ (599) NFITsII1 
9 FORMAT (215) 
DO 19 NFT=leNFIT 
READ1 (S¥TRSCA2Z11) 9 CAB (1) sCA4(1) sCA5(1) eUT9CA1(1) 
1 FORMAT (6F1005) 
READ (592) AK11sAK219AK229AK41 
2 FORMAT (4F804) 
THateut 
BK11=AK11/26 
BK22=AK22/2e 
Poesaal=1s1i 
XCAG( ITI =leShlet+AKLI*¥T(I)) 
3 TCI+l)=T¢(1)4+uUT 
Cx*# SET INITIAL VALUES 
G 
DO 8 -JvelelIl 
121 
7 W(1)=(=AK11*¥CA1L( 1) ##20) *UT 


E(1)=UT*(BKL1¥CA1 (1) **200=(AK21+AK22)*CA2(1)) 
B(1)=(AK21*CAZ(1))*UT 
CLL) =(BK22*CA2 (1) -AK41*CA4 (1) **005) ¥UT 
D(1)=(2*AK41*CA4 (1) ##0 05) ¥UT 
IF (1-3) 495915 

15 IF (1-4) 49696 

4 CAL(I+1)=CAL(1)+005*W(I 
CAZ(1+1)5CA2(1)+005*E(1 
CAB(1+1)=CA3(1)+005*B(1 
CAG(1+1)=CA4(1)+005*C(1 
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APPENDIX 4 = 
C - RUNGE = KUTTA eee (CONT'D) 
CA5(T+1)=CA5(1)+005*D(T) 
T=I+1 3 
GO To 7 


CAL(I+1)=CAL( 1) +Wwl1) 

CA2(1+1)=CA2(1)+E(1) 

CA3(14+1)=CA3(1)+B(1) 

CA4(1T4+1)=CA4(1)4+C(1) 

CA5(141)5CA5(1)+D0(1) 

1h te) 

GOMAGO! 7 

XCAL( Stl) =le/60e*(W(1)+20*W(2)4+20%*W(3)+W(4) )+CA1(1) 
XCA2Z( It ll =zle/Ge*(E(1)+20*E(2)4+20*E(3)+E(4) )+CA2(1) 
XCA3(J+1)=1le/60*(B(1)+20*B(2)+20%*B(3)+B(4) )+CA3(1) 
XCA4( Stl) ele /be*(C(1L)4+20*C(2)4+20*C(3)+C(14) 1+CA4(1) 
XCA5( Stl) ele /60k(D(1)420#D(2)4+20%D(3)4+0(4) )+CA5(1) 


RESET INITIAL VALUES 


CA1(1)=XCA1(U+1) 
CA2(1)=XCA2(U4+1) 
CA3(1)=XCA3(U4+1) 
CA4(1) =XCA4(U41) 
CA5(1)=XCA5 (U4+1) 
DO 184Js=1915 
WRITE (6910) 


LO FORMAT C'L'9/S/S/S//S/A////952Xs"TABLE 1'9/950X9' INPUT 


* DATA's/950X9!—== 


[== ----') 


WRITE (6911) AK119AK215AK22 sAK41 


Ll FORMAT (' '929Xo'K11=2!' oF 4e295X9'K22=' 9F4e298X9'K23=! 


*9F4e298X9'K44= 
Li osF4e29/929Xe'KKIL=AK1I1/2' 919X9'KK23=K23/2') 
WRITE (6914) 


14 FORMAT (' '945Xe9'INITIAL CONDITION! 9/946X9!—--=—= -- -- 


Hommmmmn ly /y52 
LXe'TIMEF00CO (HRe) 'o1l4Xs' INTERVAL H = Oe002' o/e27Xs'Al 
% sly (MOLCE/LI 
LTRE)' 95Xe'A2 =A3B =A4 =A5 =000 (MOLE/LITRE)'!'/) 
WRITE (6916) 


16 FORMAT (' '916Xe9 "REACTION! 95X9'Al'99X9!'!A2' 98X9'AZB' 9 8X 


eg 'AG' 98Xo9'AS! 
196Xe9'MATERIAL#! 92X9'TRUE Al*¥*'o9/916Xe9'TIME (HRe) 'sl6X 
% y ' COMPONENT ( 
IMOLE/LITRE) '915X9'BALNCE' 93X9'MOLE/LITRE'!) 
DO 13 125009119250 
J=I+1 
XM (J) =XCAL( IU) +2*XCA2 (9) #XCAB (I) #244*#XCA4G (4) F2*XCAS5 (4) 
T(I1)=T(1)+020007 


13 WRITE (6912) T(1) 9XCA1(U) 9XCA2 (4) 9XCAZB (4) 9XCA4( 4) 
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APPENDIX 4 = 10 


c RUNGE = KUTTA eee(CONT!D) 
*—9XCA5(4) »XM(U) XC 
LA6(1) 
12 FORMAT (! '915X9F60295X95(F8e692X) sl1X9F8eb93X9F8 0G) 
WRITE (6917) 
17 FORMAT(' 's18X9'#TOTAL NUMBER OF MOLE/LITRE EQUIVALENT 


# TO ALO PRES 
LENT IN THE REACTION SYSTEM'9/919Xs'AT A COMMON 
* INSTANEOUS TIME = A 
11 + 2*A2 + 2%A3 + 4*#A4G + 2*A5'9/918Xo' **ANALYTICAL 
* EXACT INTEGRATI 
1ON!) 
18 CONTINUE 
19 CONTINUE 
ShOP 
END 
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C PLOTTER 
fa PLOTTER 
G HK FH FE FE FE HE FE FE FE FE 36 IE 36 FE IE FE Fe HE HE IE HE HE 4b OSG HE HE 
G % * 
s * XC APLOFTERs -EXAMPLEE 2 * 
GS % % 
Cc HH HE HH HEHE OE FE I HH I He OE He EE EE FE HE HE 


C#¥*¥* THIS PROGRAM WAS WRITTEN BASED ON IBM 1627 PLOTTER 
C#e*e MANUAL 
G 


DIMENSION B(5)9C(5) 9F (5) sW(5) sCA1(5) 9CA2(5) 9CA3 (5) 
#9CA4(5) 9CA5(5) oX 
LCASITOS ) sXCAZ(105:) 9XCA34(105) sXCAG(105) oxCA5(105) 
*90(105) 9E(5)9T(19) 
19A1(19) 9A2(19) 9A3 (19) 9A4(19) 9A5(19) 
READ (592) CAL(1) 9CA2(1) sCA3(1) 9CA4(1) 9CA5(1) 
READ (592) AK119AK229AK23 9AK44 sAN11 9AN22 9AN2Z3 9ANGS 
FORMAT (8F8s6) 
NPT=101 
UT=Oe0el 
D(1)=UT 
DO 3 IT=l»NPT 
3 Di I+1)2001)4UT 
DO 8 JElesNPT 
I=1l 
7 WIT) SCMAKLTIL CALC I) ¥*¥ANI1) *UT 
E(T)SUT#(AKLI/SANLI®CA1L( I) *¥#ANILI@AK22*CA2 (1) *¥#AN2Z2"AK23 
*#CA2 (1) #¥AN2 
pe ® 
Bi T)=UT#¥(AK22*CA2 (1) *#*AN2Z2) 
CUCL) SUT#(AK23/20%*CA2Z( 1) *#ANZ3Z@AK4G4G#CAG (1) K#ANGG) 
F(T) =UT#(2#AKG4#CA4G (1) ##ANGS) 
IF (3=3) 495915 
15 IF (1"4) 49696 
~& CAL(I41)=CA1(1)+00e5*wW(T) 
CAZ(I41)=CA2(1)+00e5#E(T) 
CAB(I41)=CA3B(1)4005*B( 1) 
I) 
I) 


Nm 


CA4( IT +1) =CA4(1)4005*C( 
CAS (141)=CA5(1) +0 05#F ( 
l=I+l 
GO TRO ve 

5 CALCI+1LECA1L(1T)+wlT) 
CAZ(IT +1) =CA2(1)+E(1) 
CA3(1+1)=CA3(1)+BC1) 
CA4(I+1)=CAG(1)+C(1) 
CASI 41) 5CA5(1)4F(T) 
T=I+1l 
GOADO. TL 
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APPENDIX 4 = 12 
c PLOTTER eoe(CONT'D) 


& XCAL(Utl)=le/be*(W(1)+20*W(2)+20e*W(3)+W(4) )+CAlL(1) 
XCA2( Stl) zle/GeX (Ell) +2e%E(2)+20%E(3)+E(4) )+CA2(1) 
XCAZB(Utl)=1le/60*(B(1)+20*B(2)+20%*B(3)+B(4) )+CA3 (1) 
XCA4 (Stl) =le/60e*(C( 1) 420% C(2)+20%C(3)4+C(4) )+CA4(1) 
XCAS(Utl)=le/6e*(F(1L)+2e#F(2)+20%*F(3)4F (4) )+CA5(1) 


Cx¥**® RESET <ENITFITAL VALUES 


CA1(1)=XCA1L(U+l) 
CA2(1)2=xXCA2(U+1) 
CA3B(1)=XCA3(U+1) 
CA4(1)=XCA4(U41) 
8B CA5(1)=XCA5(U+1) 

DO 500 125910095 
J=T+1] 

500 WRITE (69501) D( I) »XCAl(U) »XCA2Z2 (4) »XCA3B (4) sXCA4S( YL) 

*#9XCA5(U) 

501 FORMAT (' '920Xs6(F100592X)) 
SX=02e625 
SY=9e5 
PP=5e 
PX=020 
PY=O06 
CALL FPLOT (190690¢) 
CALL SCALF(SX9SY 9PXoPY) 
CALL FGRID (0906 900 sPP/50e910) 
CALL FGRID (1990 9006904e05910) 
CALL FPLOT (190050e5) 
CALL SCALF (SXsSY90090e5) 
CALL FGRID (0900 s00e5sPP/5e910) 
CALL FPLOT (1910¢80e¢5) 
CALL SCALF (SX9SY9100e9005) 
CALL FGRID (391049005900¢05910) 
CALL FPLOT®4 1900 90) 
I=1l 
J=I+1 
CALL FPLOT (29D(1) »9XCA5(U)) 

140 Y=XCA5(J) 
CALL FPLOT=409D0(1) 9Y) 
IF (TEST#=D(1I)) 14191419142 


142 [=I+1 
J=I+1 
D(I+1)20(1)4UT 
GO TO 140 


141 CONTINUE . 
CALL FPLOT (190¢90e¢) 
DOr 2a. 1F "1.19 

211. READ? 4552929" AlUT) 
DO 2 Ss! l= ie. 

213 READ (59212) A2(1) 
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APPENDIX 4 = 13 
C PLOTTER eee (CONT'D) 


DO 214 IT=1919 
2E4RREAD (55212) ABI) 
DO wZ15 bali 29 

DL een C58 2 ba Aaa 
0O 216 121919 

Sl6oR EAD (5 32) 2 ASA) 

DO-217f Lalelks 

el ore AD (5.62 12 iy Fie) 

212 FORMAT (F13s7) 

CALL FCHAR (Gels 02e06900els0el 900) 
WRITE (79200) 

200 FORMAT C'REACTION TIME (HRe)*) 

CALL FCHAR (0029702087590 0l 2001 90e) 

WRITE (79201) 

201 FORMAT (REACTION PATH CALCULATED FROM KINETIC MODEL 
«* OBTAINED FROM 

1 THIS WORK") 

CALL FCHAR (2es*OellsOvels0els0e) 

WRITE “(te 202) 

202 FORMAT CYHYPOTHETICAL REACTION SYSTEMs EXAMPLE NOe 1!) 
CALL FCHAR (#103 900el75s00els0elsle57) 
WRITE (79203) 

203 FORMAT ('CONCe (MOLE/LITRE)') 

CALL FCHAR (1429003990008 90008 906) 
WRITE (79204) 

204 FORMAT ('SIMULATED EXPERIMENTAL POINTS!) 
CALL FCHAR (6282004690008 90008906) 
WRIFE (79205) 

205 FORMAT (' COMPONENT!) 

CALL FPLOT (=296090041) 
CALL POINT (2)) 
CALL FCHAR (7029004190008 90008 906) 
WRITE (79207) 
207 FORMAT. ChAZ!) 
CALL FPLOT (=296690e39) 
CALE POINT | C3) 
CALL FCHAR (729003990008 90008 90e) 
WRITE (73208) 
208 FORMAT ('A3!) 
CALL FPLOT (=296090037) 
CALL POINT (4) 
CALL FCHAR (7029003790008 90008 90e) 
WRITE (79209) 
209 FORMAT ('AG!) 
CALL FPLOT (2960290035) 
CALL POINT (5) 
CALL FCHAR (7029003590408 90008 906) 
WRITE- (79210) 
210 FORMAT ('A5*) 
CALL FPLOT (1906906) 
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€ PLOTTER gee (CONT'D) 


GAAS CER VT. 
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APPENDIX: 4 = 15 
€LEASTUSO 
AST SQ 
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% 
eo: * 
* LINEAR LEAST SQUARES POLYNOMIAL APPPROXIMATION * 
# * 
# # 
eRe ee ee ee ee ee eee 


N=NO«e OF INPUT DATA POINTS 
M=POWER® OF THE POLYNOMIAL TO BE. FITTED 


DIMENSTON. T1450) sY(50) »sSNAM(5) 

DATA SNAM/tAO =! 9tAl s=to'A2 ='94AB ='o'AG =!/ 
N219 

DO 1 Te#loeN 

READS 2 Ue ¥x( 1") 

CONTINUE 

DO 3. T=loeN 

READ S920? T4194 

CONTINUE 

FORMAT AP iS 7) 

DO” 45 1 =h'53 

Ve[+] 

VM=MéM+4] 

KKKeM+2 

WRIME “Vos KKK 

FORMAT “W* 259 /77/77//926X9'TABLE' sI2) 

WRITE (C6's69"M 

FORMAT ('O'915X9'LEAST SQUARES! 9I291Xs ORDER 
POLYNOMIALS!" ) 

CALL LEAST (TsYsNoMoaSNAM) 

CONTINUE 

tee Xl i 

END 
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SUBROUTINE REGEN 


SUBROUTINE REGEN(XsYsZsMMoN 
DEMENS LON. xX US 0) (50) 92120) 


APPEND A "4. S016 


) 


THIS SUBROUTINE REGENERATES DATA POINTS FROM THE 


RESULTED POLYN. 


WRITE VG.g:l0) 
T PGR ATt 2 4329Xs *REGENERATED 
49 Y OBSERVED" 
Pearse CY RCALCULATED Ms 3X 9! PCT 
VAR=Q6« 

HIT =Q6« 

Ni 

K2=N#2 

DO 2 T#loeN 

CAL=0 6 

DO 3. J=leM™ 
CAEaCART ACO eK 1) ee = 1) 
CONTINUE 

See bp mic NL. 
PCE=CAT/ABS(Y(1I))#1006 
VARSVAR+CAT ##2 

IF (HIMARS(PCE) 149402 
4 HT SARS (PCE? 


WwW 


DATA'//10X9'X 


ERROR! 9/) 


CoV IIE tG 9s) “Ais 1 CE PSCAls PCE 


5 FORMAT ( 9X94(F110694X)) 
VAR=VAR/ (N=MM) 
DEVEVAR ¥¥0 05 
WRITE (696) VARSDEVsHI 

6 FORMAT(//910X'VARTANCE 

%* "STANDARD DEVIATION 
1 ='9F1006//10Xs "MAXIMUM PCT 
RETURN 


iS }i2 


MEASURED’ 95x 


="59F100e¢6//10X 


ERROR =ulng eo Neo) 
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€ LEAST SQ 


G LEAST SQ 
SUBROUTINE LEAST (XeYsNsMsSNAM) 
DINENSION X50) 9Y(50) 9A(5095) sP( 20920) 9V(20) 52(5) 
*9SNAM(5) 
NPLT=0 
DO 21 IT=195 
fe Oe: Q 
MM=aiv+] 
DO 4 T=losN 
DO 4 J=leMM 
G4 ACTeJIHXCT) #¥(JH1) 
DO 5 JT=lsMM 
DO 5 J=loMM 
P(leJ)=Ue 
DO 5 K=leN 
5 P(IsV)=P(19uU)+A(KeI)#A( Kou) 
DO 6 IT=1l 99MM 
V(JT)=06 
DO 6 JY=loN 
6 VOL%EVCERY C3) XA Js) 
CALL GAUSS(PsV9MMeZ) 
WRITE (698) 
CerORMAI CUO Cte luns ite AGOEREPGIENTS OF? THE FITTED 
eR UNCTLON ARE /.) 
DO 15 T=19é™M 
Lo WRITE (6s 2) SNAMCI) sZ¢ 7) 
7 FORMAT (15X9A49F11.8) 
CALL REGEN (X9¥sZoMMsN) 
IF(NPLT) 999914 
14 CALL POLYT (X9Z9N9MM) 
9 CONTINUE 
RETURN 
ENO 
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SUBROUTINE POLYT 


SUBROUTINE POLYT(Xs29NoMM) 

DIMENSION X(50)2Z2(20) 

WRITE (691) 

POR Tre? Poe N SP LOM Teal DATAY/S/259X"%, CACCOCATED" 94x 


Be bY CALCULATE 
rots 


XMAX=O 6 

XMIN=999996 

DOr 2 T=leN 

IF (XMAX=X(1)) 39394 
XMAX=X( 1) 

TFUXRCT P=XMINDA 5°95 9 2 
XMIN=X(1) 

CONTINUE 
DELX=(XMAX=XMINI)/206 
XY=EXMIN 

DO 6 T=1920 

CAL=O06 

DO15 Yyelem™M 
CALECAL+Z (0) *XY** (U1) 
WRITE(697) XYsCAL 
FORMAT(24Xs2(F11+¢694X) ) 
KY=KY#DELX 

RETURN 

END 
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SUBROUTINE GAUSS 


SUBROUTINE GAUSS (AsRoNoX) 
DIMENSION A(20920)9R(20) 9X(20) 


C#*#* THIS SUBROUTINE CALCULATES PARAMETERS OF THE POLYNe 


M=N-1 
DO ll J=EloemM 
S=0e6 
DO 12 TUN 
U= ABS(A(I9J)) 
IF (U#S) 1291238112 
112 S=U 
Per CONTINGE 
IF (Led) 119919119 
119 DO 14 JT=uUeN 
S=A(LolI) 
A(LeT)=AlusT) 
4 AllsI)=S 
S=R(L) 
R(LIFR I) 
R(Uu)=S 
POPE ABSt(AtJIs)) }eleE=30) 1159115945 
LS WREFE (6.43) 
3 FORMAT (1H s* MATRIX SINGULAR! ) 
RETURN 
15 MM=J+1 
DO ll I=MMoN 
TFA CABS ACI ss) )—1eE=-30) lLlelllslll 
21) S=AC Qs) /A (Lod) 
AlCITsu)=CeC 
DO 16 K=MM9N 
16 AllsK)I=AlUsK)—~S*A(I 9K) 
R(T)=R(U)—S#R(T) 
11 CONTINUE 
DO 17 K=leN 
I=N+1-K 
S=0e0 
bee otk dal Ske Sd ke 
Lig MM=ay+1 
DO 18 JEMMsN 
18 S=S+A(I9JU)#X(J) 
le Se elm S LSA CT 9 1) 
RETURN 
END 
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